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ABSTRACT 

The  restoration  of  incoherent  optical  objects  which  have  been 
diffracted  by  an  optical  system  and  corrupted  by  detector  and  additive 
background  noise  is  considered.  The  approach  to  the  problem  is 
basically  numerical  and  considers  operating  directly  on  the  noisy  image 
and  point  spread  function  rather  than  the  Fourier  transform  of  these 
quantities.  The  effects  of  noise  and  the  se  of  a  priori  information 
in  the  restoration  process  are  given  particular  attention. 

Several  "optimum"  estimates  of  the  object  intensity  distribution 
are  considered.  Based  on  statistics  which  have  been  verified  in  practice, 
the  Baye's,  maximum  a  posteriori,  maximum  likelihood  and  mean  square 
error  estimates  of  the  object  intensity  distribution  are  obtained.  These 
statistical  estimates  are  compared  mathematically  and  in  many  cases 
numerically  to  other  non-etatiatical  estimates  formulated  from  control 
theory  and  dynamic  programming.  Extensive  numerical  results  have 
been  obtained  for  the  restoration  of  various  one -dimensional  objects  in 
the  presence  of  noise.  Two  monochromatic  "point  sources"  in  the 
presence  of  noise  are  shown  to  be  resolved  when  separated  by  1/5  of  the 
Rayleigh  criterion  distance.  Numerical  results  are  also  shown  for  the 
mean  square  error  as  a  function  of  a  priori  information,  the  measuring 
scheme  chosen  and  diffraction. 
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INTRODUCTION 


Statement  of  the  Problem 

All  object*  which  are  input*  to  a  physically  realisable  optical 
imaging  system  emerge  from  the  system  as  imperfect  representations 
of  the  original  object.  This  object  corruption  may  be  due  to  several 
effects  inherent  in  the  optical  system,  such  as  diffraction  and  aberration, 
as  well  as  background  and  detector  noise.  7n  fact,  the  effects  of 
diffraction  alone  can  severely  limit  system  resolution  and  effectively 
destroy,  insofar  as  the  observer  is  concerned,  die  original  object  detail. 
These  effects  are  so  pronounced  that  in  the  past,  as  Toraldo  di  Francis 
(1952)  points  out,  the  classical  Rayleigh  resolution  limit  has  been 
accepted  as  a  theoretical  limit.  However,  in  1955  Toraldo  suggested  that 
a  priori  information  about  the  object  could  be  used  to  alleviate  the 
apparent  theoretical  limits  imposed  by  diffraction. 

In  1964  J.  L.  Harris  showed,  by  using  some  well-known  results, 
that  in  the  absence  of  noise  a  priori  knowledge  that  the  object  is  of  finite 
extent  is  sufficient  to  enable  the  exact  restoration  of  a  diffracted  optical 
object.  This  work  theoretically  establishes  that  object  restoration  is 
limited  only  by  system  noise. 


l  2 

With  this  background  we  state  that  the  problem  treated  in  this  paper 
is  that  of  restoring  an  optical  object  which  has  been  diffracted  and 
corrupted  by  background  and  detector  noise.  Particular  emphasis  is 
placed  on  the  trace -off  between  restoration  detail  and  the  limiting  noises. 

Approaches  to  che  Problem 

Although  a  unique  theoretical  solution  has  been  shown  to  exist,  the 
application  of  a  working  restoration  procedure  which  performs  well  for 
a  variety  of  objects  and  optical  configurations  is  quite  another  matter. 
Basically,  three  approaches  have  been  considered  in  solving  this 
problem. 

One  approach  has  been  to  utilise  the  Fourier  transform  relationships 
that  exist  between  the  object  and  image.  Harris  (1966).  Osterberg  (1966) 
and  v/olter  (1961)  have  considered  this  approach.  However,  this 
technique  appears  to  be  difficult  to  analyse  when  applied  to  the  noisy 
observed  quantities,  and  as  yet  no  results  appear  to  be  available  which 
use  Fourier  techniques  in  such  an  environment. 

A  second  approach,  used  by  Barnes  (1966a),  is  one  of  analytically 
solving  for  the  object  by  using  the  eigenfunctions  of  am  integral  operator 
involving  an  analytical  point  spread  function.  The  drawbacks  of  this 
approach,  from  a  practical  point  of  view,  are:  (1)  the  extreme  difficulty, 
at  present,  in  staining  eigenfunctions  for  a  general  point  spread  function 
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(one  that  i»  measured),  and  (2)  the  fact  that,  to  date,  system  noise  has 
not  been  included  in  the  analytical  restoration  process. 

The  third  approach  which  is  developed  in  this  paper  is  that  of  using 
numeric  .  methods  to  perform  the  restoration.  We  reiterate  the 
practical  advantages  of  this  procedure  over  those  above.  First,  in  an 
actual  real-world  situation  only  discrete  quantities  are  available  to 
operate  on  and  effect  a  solution;  i.  e.,  the  point  spread  function  end 
image  are  usually  not  analytically  known  and  so  must  be  discretely 
measured  and  thei  operated  on  to  perform  the  restoration.  Secoi  t, 
the  observed  quantities  are  not  deterministic  but  must  be  tr»a  *>d  as 
random  variables.  This  randomness  may  be  due  to  additive  background 
fluctuations,  film  granularity,  photomultiplier  multiplicative  roise  or 
some  other  detector  noise.  These  effects  can  be  taken  into  account  by 
using  numerical  techniques  and  appealing  to  the  statistically  optimum 
minimum  mean  square  error  (MSE)  estimation  procedure. 

Flan  of  the  Paper 

The  first  section  presents  a  review  of  the  development  of  the 
image-object  relationships  and  presents  the  integral  equation  that  is  to 
be  solved  to  effect  the  restoration.  Emphasis  is  placed  on  reviewing 
tiie  assumptions  involved  in  developing  the  imaging  equation  and  in 
discussing  the  concept  of  diffraction  and  how  it  is  to  be  varied  in  this 
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paper*  Emphasis  ie  also  placed  on  the  transfer  of  the  continuous  version 
of  die  imaging  equation  to  its  discrete  representation  so  that  numerical 
methods  can  be  used  in  its  solution* 

As  a  preview  to  the  numerical  work  that  will  follow  later,  the  second 
section  after  the  introduction  reviews  previous  numerical  work  in  the 
solution  of  the  basic  integral  equation  developed  in  the  first  section. 

Also  presented  in  this  section  is  a  dynamic  programming  solution, 
developed  by  the  author,  which  utilises  an  "area"  constraint  in  the 
solution. 

The  third  section  treats  two  analytical  approaches  which  assess  the 
difficulty  in  regaining  object  information  which  has  been  obscured  due  to 
diffraction  and  noise.  Particular  emphasis  is  placed  on  the  trade-off 
between  SNR  and  diffraction  in  obtaining  this  information. 

The  fourth  section  discusses  the  noise  models  to  be  used  and 
presents  die  development  of  the  minimum  MSE  object  estimate.  This 
estimate  and  the  computational  features  discussed  in  the  last  two  sections 
constitute  the  basic  restoration  procedure  as  presented  in  this  paper. 

Results  of  performing  actual  simulated  restorations  and  of  the  study 
of  errors  involved  in  using  die  MSE  estimate  are  presented  in  the  last 
two  sections,  h’ere  again,  as  in  the  section  on  analytical  results, 
particular  emphasis  has  been  placed  on  the  ejects  of  noise  and  diffraction 
in  restoring  object  detail,  hi  addition,  emphasis  has  been  placed  on 
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presenting  actual  simulations  for  various  nolle  levels  so  that  one  may 
visually  judge  the  restoration  improvement*  These  simulations  augment 
the  important  but  somewhat  abstract  results  obtained  in  the  analytical 
section  and  serve  to  give  one  a  more  complete  picture  of  the  restoration 
procedure* 


6 


THE  IMAGING  EQUATION 

In  this  section  we  derive  end  discuss  the  relationship  between  the 
radiation  from  an  incoherent  source  intensity  distribution  and  the  corres¬ 
ponding  image  intensity  distribution.  Near  the  end  of  the  section  the 
imaging  equation  in  both  continuous  and  discrete  form  is  stated.  In 
order  to  clarify  what  is  meant  by  diffraction  as  it  applies  to  the  restoration 
problem,  the  concept  of  diffs  action  and  its  relationship  to  the  imaging 
equation  will  be  discussed.  With  this  background  the  reader  should 
be  able  to  more  clearly  assess  what  is  meant  by  object  restoration  in 
the  presence  of  diffraction. 

The  Optical  Configuration 

The  optical  configuration  to  be  used  in  the  derivation  is  shown  in 
Figure  1.  However,  it  should  be  mentioned  that  the  imaging  equation 
developed  here  with  the  subsidiary  definitions  necessary  to  utilise  this 
equation  can  be  derived  from  several  optical  configurations.  Two  such 
derivations  are  mentioned  below. 

Kelstrom  (1964)  derives  this  equation  with  suitable  geometric 
assumptions  which  make  it  possible  for  a  lens  or  absorber  to  replace 
the  aperture.  His  approach  is  quite  natural  if  one  wants  to  develop  the 
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practical  aspect  of  an  instrument  which  observes  objects  at  great 
distances  from  the  aperture* 

Cutrona  and  others  (I960)  use  a  configuration  composed  of  appropri¬ 
ately  spaced  lenses  which  is  often  used  in  the  laboratory.  Using  this 
system,  the  only  assumptions  necessary  to  derive  the  imaging  equation 
are  that  the  lenses  together  serve  to  image  the  object  plane  onto  the 
image  plane  and  that  the  lenses  are  aberration  free. 

The  Scalar  Theory 

To  begin  with,  the  disturbance  in  the  object  plane  is  taken  as  an 
electromagnetic  field  at  optical  frequencies.  The  complete  derivation 
of  the  imaging  equations  by  finding  appropriate  solutions  of  Maxwell's 
equations  has  been  accomplished  only  in  a  few  idealised  situations 
(Stone,  1963).  However,  an  approximate  theory  has  been  developed 
which  allows  one  to  solve  the  imaging  relations  for  a  variety  of  aperture 
shapes.  The  approximate  theory  treats  the  electromagnetic  disturbance 
as  a  scalar  field  and  neglects  effects  due  to  polarisation.  Fortunately, 
the  imaging  equations  developed  from  this  theory  predict  the  associated 
distributions  to  high  accuracy  when  certain  restrictive  although  not 
impractical  assumptions  are  made. 

The  wave  equation,  which  the  amplitude  of  the  wave  in  the  scalar 


theory  must  obey,  is 
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2  .  i  a  ii  „ 

v u+  T  2  ■ 0 
c  a  t 


where  u  is  the  scaler  amplitude  of  the  wave  and  e  is  the  free  space 
velocity  of  light.  The  monochromatic  (wavelength  x)  solution  of  this 
equation  in  spherical  coordinates  may  be  written 


u(R,  t)  ■  ““»  cos  [k(R-ct)  +  0] 


where  a  is  the  constant  fixing  the  amplitude,  R  is  the  spherical  radial 

i 

variable,  k  *  2 w/x»  and  0  is  a  phase  angle.  This  solution  represents 
the  wave  propagated  from  a  point  in  the  object  plane. 

It  is  desired  to  consider  many  such  points  located  at  (or,  (3 )  in  the 
object  plane,  and  so  we  write  for  the  object  amplitude  distribution 


u(R,a,p,t)  ■  c  R|  [A(a,?,t)e  ] 


where  A(ar,p,t)  is  the  complex  amplitude,  R  t means  "real  part  of,  " 
and  c  is  a  constant  which  takes  the  inverse  distance  relationship  into 
account. 


. . I . . 
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Here  we  make  the  assumption  as  discussed  by  Helstrom  (1964)  that 
the  observation  time  is  to  be  small  compared  with  the  reciprocal  band¬ 
width  of  the  radiation  and  that  no  information  regarding  the  instantaneous 
time  dependence  will  be  used*  Dropping  the  time  dependence  and  for  our 
purposes  the  unimportant  constant  c,  we  have 

«(*.«.*)  ■Afa.Pje11*.  (4) 

Use  of  Huygen' s  Principle  and  the  Geometric  Assumptions 

i 

4? 

The  Fowler  transform  relations  which  will  be  developed  are 
essentially  those  established  by  Huygen,  Fresnel*  and  Kirchoff.  For 
a  more  detailed  treatment,  the  reader  is  referred  to  Stone  (1963)*  O'Neill 
(1963),  Bora  and  ’Jolf  (1959),  and  Be  ran  and  Parrent  (1964). 

Consider  ;  wave  due  to  a  point  source  in  the  object  plane  which 
propagates  toward  the  aperture.  In  essence*  Huygen1  s  principle  states 
that  once  this  wave  has  propagated  to  points  which  are  just  through  the 
aperture  opening,  then  these  points  themselves  may  be  regarded  as 
sources  of  secondary  wavelets  that  again  propagate  with  the  usual  solution 
(2)  for  a  point  source.  The  disturbances  due  to  all  of  these  "secondary" 
sources  which  reach  the  image  plane  then  may  be  appropriately  summed 
to  obtain  the  image  plane  distribution.  Finally*  to  obtain  the  image 

* 

distribution  due  to  several  point  sources  or  a  continuum  of  point  sources 

i 
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in  the  object  plane  we  merely  amn  or  integrate  over  the  object  plane 
coordinate*. 

Following  the  above  procedure  we  seek  to  find  the  amplitude  at  an 
arbitrary  point  juit  to  the  right  of  the  aperture  opening.  At  a  point 
to  the  left  of  the  aperture 

utRj.e.pj-AKJlje1**!.  (5) 


Dealing  only  with  the  exponential, 

IkR,  i*R10  lk(R,-R10) 

e  a  e  e  . 


<*> 


Referring  to  Figure  1  the  length  R  from  Oj  to  O  can  be  written  in  two 
uaya: 


R  aRj2-  fr-e)2  •  («T^)2  a  RJ()2  -«2-p2, 


(7) 


alao 


Rj2  -  R102  -  (Rj-R^HRj  +  R10)  «Y2  +«r2-2(fl7  +<rfJ).  (8) 


Thua 


R  u  _  yi  '♦‘y2  2fya  +  yp ) 

1  10  R,  +  R10*  R,+Rl0  • 


We  now  disregard  the  squared  term  in  the  exponential  on  the  grounds  that 
the  Fraunhofer  far -field  conditions  hold  and  that  the  open  portion  of  the 
aperture  plane  is  small.  However,  we  still  make  the  assumption  that 
tiie  dimensions  of  the  aperture  are  many  wavelengths,  an  assumption 
which  is  made  in  order  that  the  scalar  wave  theory  will  agree  with 
experimental  results  (Stone,  1963),  We  also  assume  that  the  source  is 
sufficiently  Ur  removed  from  the  apertiure  so  that  Rj  +  R1(J  «  2R  in  the 
denominator.  Thus,  the  electric  field  just  to  the  left  of  the  aperture  is 


tkR,„  -ikCliHI 
1#  •  R  . 


In  order  to  simplify  the  above  expression  (10)  the  following  spatial 
frequency  components  are  defined: 


f  «  -JL 

y  *R  • 


f  ■  -~ 
«■ 


Expression  (10)  becomes 
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To  obtain  the  total  electric  field  due  to  a  point  source  in  the  object 
plana  we  must  integrate  over  all  of  the  ’'secondary"  point  sources  in  the 
aperture  plane*  Thus,  the  electric  field  in  die  image  plane  due  to  a  point 
source  in  the  object  plane  is 


Thus  far  we  have  dealt  entirely  with  the  amplitude  the  electro¬ 
magnetic  field*  However,  throughout  the  derivation  that  follows  it  will 
be  assumed  that  we  are  dealing  with  incoherent  light  and  that  the  measured 
quantities  are  intensities,  which  add  linearly*  Coherent  light  could  also 
be  assumed  throughout  and  the  form  of  the  imaging  equation  would  be  the 
same  as  that  derived  for  the  incoherent  case*  In  this  respect  the  numerical 
techniques  used  could  be  applied  to  e idler  coherent  or  incoherent 
illumination,  but  the  simulated  results  in  the  last  section  apply  only  to 


the  incoherent  case* 
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The  intensity  (watts)  at  {£,  r\)  due  to  a  point  source  at  (a,p)  is 

*  I  l2 

■ya  0(6*n)y  •  ly^te.^)  I  (*?) 

where  the  star  as  used  here  denotes  complex  conjugate.  Performing 
this  multiplication 


b<a,IU*Tl)«  ffff  U(*,|3)  I  G fl'.OQdJ,/) 

W' 

-i2ir  (0(^-4)+ 


•  6 


+i2ir(4(fY-fy)+nifr-fr)3 


•  e 


i  i 


dfydfydf^d f,.  (18) 


To  obtain  the  total  intensity  due  to  a  continuum  of  point  sources 
in  the  object  plane  we  integrate  over  the  object  plane*  Thus 


00 

b(£,Ti>  -  JjMo.fi.i.n)  <todP, 


(19) 


and  using  (18) 


rrm  *  .  .  + 

M6.  ^1)  •JJ  flCfy.  V  G  «y,  y  e 


r  rn  12  +p^-<)]  -1  ,  . 

*  /  /  |A(nr,.8)J  e  dfydfydf^df^. 


(20) 
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The  integral  in  bracket*  is  defined  as  the  spatial  frequency  spectrum  of 
the  object  intensity  distribution*  The  quantity  1a(o#!3)|2  is  defined  as 
the  object  intensity  distribution  and  will  be  denoted  x{c  4  3 ).  The  integral 
in  brackets  is  also  the  two-dimensional  Fourier  transform  of  x(cr,3 ).  We 

t  t 

denote  it  as  ^(fy-fy,  f^-f^).  hi  general,,  capital  letters  will  denote  the 
'  patial  frequency  transform  of  the  analogous  function  in  lower  case  letters. 

i  t 

To  facilitate  further  simplification  let  ty  a  fy-fy  and  t^  a  f^-f^, 

it  i  i 

Thus  £y  *  £y-tyi  dfy  a  -dty,  and  u  f  -t  ,  d£ff  a  -dt^,  and 

»  i2*[£ty  +  Tit J 

*  If  X(ty,tr)  e 


(21) 


•  /foty.yaVv  W***, 

-^co 


dfcydt^. 


The  integral  in  brackets  above  is  defined  as  the  Fourier  transform 
of  the  point  spread  function.  It  is  to  be  denoted  as  K(ty,  t^ ).  It  will  be 
assumed  throughout  for  the  configuration  considered  that  this  same  point 
spread  function  applies  for  a.'l  point  sources  in  the  object  plane.  This  is 
the  spatial  invariance  assumption  one  uses  in  an  "idealised"  optical 
configuration  (O* Neill,  ly63).  With  this  definition 


b(5.r,)  -  /J’x(ty,tir)H(tv,t(f)« 


i2ff  [§ty  +  rjtjJ 


dtydt.  (22) 


00 


47 


Taking  the  Fourier  transform  of  both  side  of  the  above  equation, 


)-H(fy,f  )X(fy,f  ). 


(23) 


By  the  convolution  theorem 


60 

*  ffm^.  q-0)  3c(afp)  da  dp, 


•v. 


The  one-dimensional  version  of  (24)  is 


(24) 


M4)  *  J x(o)da, 

-00 


(25) 


Since  it  is  assumed  that  the  object  is  of  finite  spatial  extent.  Equations 
(24)  and  (25)  may  be  rewritten  as 


Hi.  T|)  -  If  M£->,  )  x(a,  P  )<todP, 

R2 

a 

b(4)  -  [Mi**)  x(«)d*. 


(26) 


(27) 


where  R,,  denotes  the  region  in  two-dimensional  object  space  where 
x(«,P)  is  nonzero. 
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It  if  our  objective  to  solve  the  integral  Equations  (26)  and  (27) 
for  the  object  when  the  image  and  point  spread  function  are  specified* 
These  equations  represent  the  diffraction  of  an  object  by  an  optical 
system  and  are  known  as  the  one  and  two  dimensional  versions  of  the 
imaging  equation*  Of  course,  it  is  obvious  to  tibia  point  that  no  noise 
is  represented  in  Equations  (26)  and  (27)*  Noise  effects  are  treated  in 
a  later  section* 


The  Concept  of  Diffraction 

when  the  point  spread  function  is  much  larger  in  spatial  extent 
than  the  object,  then  the  image  is  nearly  equal  to  the  point  spread 
function.  This  is  a  case  of  large  diffraction,  and  the  image  bears  little 
resemblance  to  the  original  object*  Alternatively,  when  the  point  spread 
function  is  small  spatially  in  comparison  with  the  object  the  image  is 
nearly  e^oal  to  the  object  and  the  diffraction  effects  are  almost  negligible. 

In  order  to  provide  a  measure  of  diffraction  we  define 

n  ^  Point  spread  extent  -2_. 

*  Object  extent  * 

where  by  point  spread  extent  we  refer  to  the  "Airy  disk"  of  the  point 
spread  function  and  by  object  extent  we  refer  to  a  measure  of  the  outside 
size  of  the  object  (where  the  object  is  nonzero);  i.  e.,  the  size  of  a  two- 
point  scurce  object  is  taken  as  the  separation  between  the  points. 


The  ratio  R  may  be  varied  in  two  ways*  We  can  either  fix  the 
aperture  and  vary  the  object  or  fix  the  object  and  vary  the  aperture* 

Both  methoda  of  varying  R  produce  the  aa me  qualitative  results.  In 
this  paper  we  have  considered,  so  far  as  the  simulated  results  are 
concerned,  a  fixed  optical  system  which  views  v  irious  objects*  This  is 
generally  the  case  treated  in  practice  since  the  cause  of  a  "diffraction 
limited  system"  is  usually  due  to  the  fact  that  either  the  aperture  was 
fixed  when  the  measurements  were  made  or,  on  realising  that  diffraction 
was  severe,  the  aperture  was  opened  as  far  as  possible  but  diffraction 
was  still  apparent. 

Since  we  are  presenting  results  for  a  fixed  optical  system  that  views 
various  objects,  it  was  convenient  to  fix  the  aperture  in  reference  with 
the  classical  Rayleigh  criterion  for  resolving  two  point  sources.  The 
Rayleigh  criterion  resolution  distance  was  set  at  unity,  and  thus 


R  * 


2*0 

object  size 


(29) 


and  R  n  2. 0  corresponds  to  the  spacing  of  two  point  sources  which 
Rayleigh  proposed  were  resolved*  However,  as  we  shall  demonstrate, 
two  sources  can  be  resolved  even  in  the  presence  of  noise  when  R  »  10,0. 


Now  that  the  concept  of  diffraction  has  been  reviewed  and  put  in 
context  with  regard  to  the  restoration  problem,  we  consider  next  the 
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presentation  of  +He  numerical  technique  to  be  used  in  solving  the  imaging 
equation* 


The  Numerical  Technique 

We  desire  to  solve  the  imaging  equation  by  a  numerical  technique* 
Thus,  we  are  led  to  consider  approximate  discrete  representations  of  the 
various  functions  involved*  The  one -dimensional  version  (27)  for  a  single 
point  ^  may  be  written  as 


V*  V-j’Vi  ,30) 

where  b^  »  b(^i),  h^  *  h{£i#  a/),  x^  a  x(or^),  w.  is  the  quadrature  weight, 
and  «  ^  is  the  quadrature  error  involved  in  transferring  from  continuous 
to  discrete  quantities* 

The  quadrature  error  is  very  important  in  the  restoration  process, 
particularly  if  it  is  assumed  that  the  object,  point  spread  function  and 
image  are  all  continuous  functions*  The  quadrature  error  can  be 
regarded  as  depending  primarily  upon  the  number  of  points  (M)  chosen 
and  the  type  of  quadrature  weight  (w^)  used*  Evidence  indicating  the 
proper  choice  of  these  quantities  for  the  one-dimensional  case  is 
presented  in  a  later  section*  Thus,  at  present,  we  neglect  the 
quadrature  error  and  defer  this  discussion  to  a  later  section. 
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Generalizing  (30)  to  N  image  point*  £j»  i^****^'  we  **&ve 
N  equations  of  the  form  (30)  which  may  be  written  in  matrix  form  as 

b  »  Ax  (31) 


where  the  elements  a 


ij 


in  the  A 


matrix  are  defined  a* 


*  h..w. 
ij  J 


i  «  1....N 

j  ■  1,  •  •  •  M 


and  b  and  x  are,  respectively,  Nxl  and  Mxl  dimensional  column  vectors. 
The  two-dimensional  version  (26)  may  also  be  represented  in  the  form 
of  (31)  by  a  suitable  definition  of  the  points  in  two  dimensional  space. 
This  definition  is  illustrated  in  the  Appendix. 

Equation  (31)  is  the  discrete  numerical  representation  of  Equations 
(26)  and  (27)  used  throughout  the  paper.  Observation  of  (31)  in  light  of 
the  objective  (that  of  solving  for  x  given  A  and  b)f  the  obvious  solution 
for  the  case  when  M  «  N  and  A  is  nonsingular  is 

x  a  A-1b  (32) 

where  A  '  is  the  inverse  matrix  for  A. 
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However,  such  a  simple  solution  is  essentially  of  no  value  when 
applied  to  the  practical  optical  problem*  The  two  outstanding  difficulties 
experienced  by  the  author  and  others  (Phillips,  1962;  Bellman,  1965) 
are  stated  below*  First,  in  arriving  at  (32)  it  is  assumed  that  no  noise 
is  present  in  the  system*  The  use  of  this  assumption  is  very  important 
since  excessive  accuracy  may  have  to  be  used  in  order  to  obtain  the 
correct  solution*  Se .  tnd,  as  the  diffraction  becomes  appreciable,  the 
A  matrix  approaches  a  singular  matrix,  and  from  a  computational  point 
of  view  the  successful  inversion  of  the  A  matrix  is  nearly  impossible. 

Both  of  these  difficulties  can  be  alleviated  by  taking  the  noittv,  and  the 
a  priori  information  in*  account*  This  fact  is  demonstrated  in  the 
simulated  restorations  presented  in  a  later  section. 

Phillips  (1962)  has  also  developed  a  technique  based  on  control 
theory  which  overcomes  the  problems  mentioned  above*  Because  of  the 
similarity  of  his  solution,  as  developed  further  by  Twomey  (1963),  and 
the  minimum  mean  square  error  (MSE)  solution,  in  the  next  section  a 
review  of  their  approach  to  the  problem  is  presented*  Bellman  (1965) 
has  further  extended  the  computational  procedure  developed  by  both 
Phillips  and  Twomey  and  has  applied  dynamic  programming  to  the  solution 
of  the  problem.  Hence,  we  shall  also  consider  the  dynamic  programming 
approach  and  attempt  to  evaluate  its  merits  and  limitations  as  it  applies 
to  the  restoration  problem. 
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NUMERICAL  SOLUTIONS  AND  METHODS 
IN  THE  ABSENCE  OF  NOISE 

Review  of  the  Control  Theory  Approach 


Because  the  final  form  of  the  solutions  developed  by  Phillips  (1962) 


and  Twomey  (1963)  is  very  similar  mathematically  and  conceptually  to 
the  MSE  solution  derived  later,  a  rather  extensive  review  of  these  papers 
is  presented  in  this  section* 

Recall  the  one -dimensional  version  of  the  imaging  equation  (27). 
Phillips  (1962)  mentions  that  this  integral  equation  of  the  first  kind  can 
be  unstable  in  that  infinitesimal  changes  in  b(£)  can  cause  large  changes 
in  x(or)  and  that  the  success  in  solving  this  equation  by  any  method  depends 
largely  on  the  accuracy  of  b(g)  and  the  shape  of  the  kernel  h(( -*-')• 

Since  the  solution  will  depend  on  the  accuracy  of  b(£),  Phillips 
suggests  that  (27)  be  altered  to  read 


a 

b($)  +  e(|)  *  ^h(|-o)x(a)  do  (33) 

where  e(£)  is  an  unknown  arbitrary  bounded  function.  The  solution  to 
this  equation  is  not  unique,  but  now  we  seek  to  find  the  best  solution 
from  a  family  of  solutions.  The  "best  solution"  is  referred  to  in  the 


Wirt||ti||i||riiir|fili||||yilTi^||[||)i  iiiWlflSMWiiitfiinir  i— jiwwwm . ******  ■>»«*-**«**#* 
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sense  that  further  constraints  are  placed  on  the  solution  which  enable 
one  to  solve  the  problem* 

Phillips  introduced  a  smoothness  constraint,  which  is  that  the 
quantity 

00 

/[*"(«>]  ^ %  <34> 

“00 

be  small,  where  x"(£)  denotes  the  second  derivative  of  x(£).  Numerically 
(34)  may  be  approximated  by 


ski  ***1 +  *i+i)2,  (35) 

i»l 


Phillips  then  considers  the  function  e(£)  to  be  bounded  in  the  following 


way: 


(36) 


or  equivalently 


M  2 

S  e  <  E.  (37) 

iel 


9 


Now  we  can  prop,  e  the  following  minimisation: 
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M  2 

rrdn  2  (x  -2x  +  x.+1)  * 

<V  i-1  1  1  i+1 


subject  to 


(38) 


N  M  2 

2(2  a  x  -b  f  -  E  ■  0.  (39) 

ieljel‘J  3 

Follov/ing  the  general  minimisation  procedure  using  -^grange  multipliers 
the  problem  posed  in  (38)  and  (39)  is  equivalent  to  seeking  the  (  xj 
which  minimise 


M  2-1 
R(x)  a  2  (xi_1-2x.  +  xi+1)  +x 

i«l 


N  M  -  2 

-v  -E 


(40) 


where  is  the  Lagrangian  multiplier  (C  our  ant  and  Hilbert,  1953; 
Hildebrand,  1952). 

If  we  continued  on  with  the  Lagrangian  method  we  would  find  partial 
derivatives  of  R(x)  with  respect  to  the  x^  and  k"1.  From  these  equations 
one  then  solves  for  the  x^  and  k"1.  In  this  case,  however,  we  cannot 
solve  explicitly  for  k  *  unless  the  bound  E  is  specified.  Thus,  to  make 
the  problem  more  mathematically  tractable  we  consider  the  quantity 
k”1  to  be  known  and  E  to  be  unknown.  This  is  often  done  in  practice, 


i 

&!• 

VH. 


M 


p- 
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and  one  can  reason  that  the  parameter  \  is  to  be  varied  until  some 
specified  bound  is  attained  (Bellman,  1962)* 

When  1*  nown  we  can  equivalently  choose  x^  to  minimise 

N  M  M  , 

R(x)  -  E  (Ea..x  -b.f  +  k  £  (x.  j-2x.  +  x.  )  .  (41) 

i«l  J  *  isl  1-1  *  i+1 

Now  the  second  term  can  be  regarded  as  the  constraint,  while  we  seek 
to  minimise  the  squared  error. 

Note  that  if  x  is  allowed  to  be  negative  the  R(x)  can  be  made  zero 
regardless  of  (x^).  Thus  only  non-negative  values  of  ^  are  to  be 
considered. 

To  find  (xj  we  differentiate  R(x)  with  respect  to  and  equate 
the  result  to  zero.  Thus, 


mm 


N  M 


«  2 

ial 


(42) 


+  (xk-2"  ^-1  +  6V4xk+l+Xk+2>  3  ° 


where,  x  a  xw, ,  =  0,  k  »  3,4,  .  ..M-2.  In  matrix  form  this  equation 
o  M+1 

can  be  written 


A' Ax  -A'b  +  t^Hx  =  0 


(43) 
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0  0  . 
0  0  . 
1  0  . 
0  1  “4  6  -4  1  t 


The  prime  is  used  to  indicate  a  matrix  or  vector  transpose.  Thus, 
the  minimizing  vector  x  is 


urhofA 


W 


5 

•4 

1 


-4  i  0 
6  -4  I 

-4  6-4 


x  a  (A»A  +  *Hf 1  A«b.  (44) 

This  is  the  result  obtained  by  Twomey  (1963)  and  is  a  slightly  modified 

e 

form  of  the  original  expression  developed  by  Phillips  (1962). 

From  (44)  we  see  that  if  \  is  zero  the  solution  reverts  back  to  the 
overdetermined  pseudoinverse  solution  (M<N)  or  tint  inverse  solution 
(M=N).  Values  of  ^X)  weight  the  constraint  more  heavily  and  arc  chosen 
in  accordance  with  the  amount  of  "smoothing"  necessary  to  control  the 
unstability  of  the  system. 

Once  the  solution  vector  x  is  found  the  individual  "errors"  (e/j  may 
be  determined  from 


e  a  Ax  -  b. 


(45) 


These  values  may  be  used  as  a  criterion  t  ■  choose  x  (Phillips,  1962). 
For  k  «  0,  t\«s  e^  are  zero. 
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The  choice  of  a  suitable  value  for  x  may  lead  one  to  be  rather 
skeptical  about  the  use  of  such  a  parameter*  Experience  has  shown, 
however,  that  in  applying  this  technique  to  the  optics  problem  the 
ill-conditioning  of  the  original  system  can  lead  to  solutions  which  conflict 
with  the  a  priori  knowledge  that  the  incoherent  object  intensity  distribution 
cannot  be  negative  or  exhibit  large  positive  or  large  positive  and  negative 
oscillations*  Based  on  this  rather  limited  prior  information,  a  practical 
choice  of  x  can  be  made* 

In  practice  we  can,  by  obtainfng  a  solution  for  several  values  of  x, 

determine  a  minimum  value  \  .  so  that  the  solutions  do  not  conflict 

min 

radically  with  this  a  priori  knowledge*  For  a  range  of  values  above  this 
minimum,  the  solutions  do  not  change  appreciably,  and  thus  are  all 
"acceptable  solutions."  With  this  knowledge  and  the  use  of  further 
computational  schemes  we  can  from  a  practical  point  of  view  circumvent 
the  arbitrariness  in  the  choice  of  x> 

Other  Constraints 

Twomey  (1963)  and  Bellman,  et  al.  (1964,  1965)  utilize  other 
constraints  which  are  applicable  to  their  par  ticular  problems.  Twomey 
(1963)  suggested  a  constraint  which  minimizes  the  sum  of  squares  of 
the  differences  between  the  actual  solution  vector  x  and  an  a  priori  vector 
p*  In  this  case  we  seek  to  minimize 
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N  M  ?  M  2 

R(x|  +  •  <46» 

The  vector  solution  (obtained  in  the  same  straightforward  way  as  bsfore) 
is 


x  a  (A«A  +  sl)ml  (A«b  +  *p)  (47) 

where  1  is  an  MxM  identity  matrix*  Twomey  states  that  a  general  form 
of  the  solution  when  many  different  constraints  are  considered  is 

x  a  (A«A  +XH)"1  (A«b  +  m>).  (48) 

where  H  is  a  matrix  and  p  is  a  vector,  and  they  are  to  be  specified 
by  the  constraint  used. 

Here  we  call  attention  to  the  fact  that  (47)  and  (48)  are  very  similar 
mathematically  to  the  MSE  solutions  (149)  and  (155)  which  are  derived 
later*  Conceptually,  the  methods  used  in  deriving  these  equations  are 
also  very  similar  in  that  the  basic  motivation  in  each  case  was  to  utilize 
a  priori  information  in  producing  a  better  solution.  Further  analogies 
between  the  MSE  and  control  theory  approach  are  deferred  to  the  Object 
Estimation  Section. 


I 


30 


Sequential  Approximations 

The  method  of  sequential  approximations  or  iteration  shown  here  was 
first  mentioned  by  Bellman,  et  al. ,  (1965),  It  was  found  to  be  very 
useful  in  the  restoration  problem  and  was  used  extensively  in  obtaining 
the  simulated  restorations  presented  later*  This  method  was  applied  to 
both  the  MSE  solution  (15 r*)  and  the  control  theory  solution  (47).  The 
method  consists  simply  of  continuously  replacing  the  a  priori  mean  or 
vector  by  an  updated  version  as  shown  below: 

x  s  (A'A  +  XI)"1  A'b  +  (A'A  +  xl)"1  x  (49) 

p  n-i 

This  is  a  sequence  in  which  x&  approaches  the  vector 

x  =  A-1h  (50) 

for  A  nonsingular*  The  convergence  of  this  sequence  is  seen  as  follows. 
For  convenience  we  define  the  matrix 


B  **  (A'A  +  xl) 


(51) 


31 


and  the  vector 

d  s  A«b,  (52) 

Now  for  x>0  and  A  nnsingular  the  eigenvalues  of  xB  are  all  less  than 
unity  (Bellman,  et  al. ,  1965).  To  show  that  the  sequence  (49)  converges, 
we  note  that  it  is  a  Cauchy  sequence 

Xn  *  \-l  *  *B(xn-l  -  Xn-2> 

*  '  Xn-3>  (53! 

xn  "  xn-l s  Bd  +  xo(I  +  xBH  * 

whence,  as  n  approaches  infinity,  approaches  zero  and  and 

x^_j  approach  x.  Thus 

x  a  (A'A  +  xl)”1  A'b  +  x(A'A  +  xlf1*  (54) 

and 


x  *  A  *b 


(55) 


32 


Again  we  stress  that  unless  we  can  successfully  invert  A  (from  a 
computational  point  of  view)  and  possess  essentially  deterministic  know¬ 
ledge  of  the  vector  b,  then  the  vector  x  will  not  correspond  to  the  object 
vector  x»  Procedures  have  been  designed  into  the  solutions  which  alleviate 
the  lack  of  this  knowledge*  However,  this  does  not  minimize  the  importance 
of  the  iteration  technique,  since  its  use  allows  one  to  effectively  vary  the 
a  priori  information  in  a  much  easier  way  than  by  inverting  a  matrix  for 
each  new  value  of  x* 

Dynamic  Programming  Solutions 

Review  of  dynamic  programming 

To  present  some  continuity  with  regard  to  dynamic  programming 
procedures,  the  following  brief  review  is  presented.  The  reader  is 
referred  to  Bellman  (1957,  I960,  1962)  for  further  details.  A  classically 
simple  yet  exhaustive  example  involving  the  principles  of  dynamic 
programming  would  serve  this  purpose,  bat  such  an  example  is  hard  to 
find  for  the  procedure  described. 

Dynamic  programming  is  a  term  used  to  describe  the  mathematical 
theory  of  performing  a  sequence  of  decisions,  or  more  formally,  the 
theory  of  multistage  decision  processor  In  this  paper  it  is  used  as  a 
mathematical  tool  to  sequentially  compute  the  vector  components  of  the 


object  vector  x. 


33 


The  'ajor  advantage  of  dynamic  programming  over  the  solutions 
just  considered  is  that  no  matrix  inverses  are  needed.  Thus,  this 
procedure  can  alleviate  matrix  inversion  errors  for  large  or  nearly 
singular  matrices.  However,  the  procedure  also  has  limitations  for  the 
problem  being  considered.  One  such  limitation  is  due  to  the  use  of  time 
as  an  artificial  index  of  available  information,  a  facet  to  be  discussed  latev 
in  this  section. 

Dynamic  programming  solutions  inherently  depend  upon  the  optimality 
principle,  which  is:  An  optimal  policy  has  the  property  that  whatever  the 
initial  state  and  initial  decision  are  the  remaining  decisions  must 
constitute  an  optimal  policy  with  regard  to  the  state  resulting  from  the 
first  decision.  The  general  proof  of  this  principle  is  proved  by  a 
combination  of  induction  and  contradiction.  A  specific  proof  is  found  in 
Bellman  (1962). 

For  a  more  detailed  review  we  consid'  u*» crete  deterministic 

process;  deterministic  in  the  sense  that  the  result  of  a  decision  is  uniquely 
determined  by  the  decision  and  discrete  in  the  sense  that  the  process 
consists  of  a  finite  numbe.  of  stages. 

Bellman  (1957)  defines  a  state  vector  p  s  (p^,  wkid*  is  a 

member  of  a  set  D  and  a  sequence  of  transformations  T  *  {T^}  .  The 
transformations  have  the  property  that  p  c  D  implies  Tq(p)  t  D  also. 

That  is,  the  transformations  have  the  property  of  transforming  the  state 
vector  at  any  stage  of  the  process  into  its  original  set. 
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A  policy  is  any  choice  of  the  set  of  variables  which  yields  an 
allowable  sequence  of  decisions.  An  optimal  policy  is  the  choice  of  q 
q^»  which  in  our  case  minimizes  a  preassigned  function  of  the 
final  state  p^«  Corresponding  to  each  we  have  T  »  and  thus  we  may 
equivalently  regard  a  policy  as  a  selection  of  transformations.  The 
preassigned  function  R(p^)  of  the  final  state  is  denoted  as  the  criterion 
function. 

Now  we  define 


*N(P)  a  niin  R  <pN)  (56) 

where  f^(p)  is  the  N  stage  return  obtained  starting  from  an  initial  state 
p  and  using  an  optimal  policy. 

We  can  use  the  optimality  principle  to  obtain  a  recurrence  relation¬ 
ship  for  the  functions  of  the  set  f^(p)  *  Suppose  a  transformation 
is  chosen  as  a  result  of  the  first  decision.  The  new  state  vector  is  T  (p). 

q 

The  minimum  value  of  the  criterion  function  as  a  result  of  the  next  N-l 
stages  is  f^_^(T^(p))  from  (56)  above.  Using  the  optimality  principle,  if 
it  is  desired  to  minimize  the  total  N-stage  criterion  function,  q  must  be 
chosen  so  that 


tN(p)  =  min  fN_,  (Tq(p» 


(57) 
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for  N  >  2,  with 
•-*  9 


fj(p)  S  min  Tq(p).  (58) 

Equations  (57)  and  (58)  ara  the  crucial  relationships  used  in  discrete 
deterministic  dynamic  programming.  Once  these  relationships  are 
established  for  the  particular  problem  at  hand,  the  solution  is  well  on  its 
way  to  completion. 

Stochastic  dynamic  programming 

A  similar  formulation  is  available  for  a  discrete  stochastic  process 
in  which  the  basic  difference  is  that  a  decision  results  in  a  distribution 
of  transformations.  Several  additive  noise  forms  were  considered  under 
this  formulation,  but  for  the  forms  considered  only  first  moments  of  the 
additive  noise  entered  directly  into  the  solution,  while  second  moments 
and  first  moments  entered  into  the  final  error  function.  Since  we  can 
assume  without  loss  of  generality  that  the  noise  mean  is  zero,  these 

solutions,  except  for  the  error  term,  were  the  same  as  the  deterministic 
solutions. 

Deterministic  solution  with  area  constraint 

Bellman,  et  al.  (1964,  1965)  has  obtained  dynamic  programming 
solutions  for  two  constraints,  one  which  involves  the  use  of  an  a  priori 
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vector  p  as  in  (47)  and  another  wKch  utilizes  a  smoothness  constraint. 
Here  we  obtain  a  solution  which  uses  an  area  constraint  and  is  suited  to 
the  a  priori  information  available  in  the  problem. 

We  first  show  that  the  area  under  the  object  is  attainable  from  the 
image  and  point  spread  function.  The  area  under  the  object  distribution 
is  the  product  of  the  areas  under  the  image  and  point  spread  functions. 
This  is  easily  demonstrable.  Consider  the  one-dimensional  version  of  the 
imaging  Equation  (25)  which  is 

b(|)  *  h(g)*x(£)  (59) 

where  the  star  as  used  here  denotes  convolution.  The  analogous 
exprer  jion  in  the  spatial  frequency  domain  is 

B(f)  »  H(f)X(f)  (60) 


where 


Expression  (60)  holds  for  all  f,  thus  when  £  =  0,  we  have 


00  00  00 
J b(uas  =  j*K  )<H  J  h(|)d| 

"co  “oo  "00 


(61) 
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which  is  the  desired  product  of  areas.  The  numerical  approximation  for 
these  integrals  yields 

M  N  N 

c  =  =  x  w  a  2  b  s  /  2  h  v  (62) 

i=l  11  ial  1  1  U1  1  1 

where  w^,  s^  and  are  quadrature  weights. 

Since  c  is  attainable  from  (62)  we  can  use  (62)  as  the  constraint  and 
minimize 


V*! 


M  2 

a  *(  £  W  X  -c) 

i=l 


N  M  2 

+  2(2  a.  .x.-b.)  . 
,  »JJ  1 


(63) 


Before  proceeding  with  the  dynamic  programming  solution,  we  shall 
consider  the  vector  solution.  Differentiating  R^fx)  with  respect  to  x^ 
and  equating  the  result  to  zero. 


aVx) 


M 

.  MS  Wft-C)  wk 
ial 


N 
+  2 
ial 


M 

a..  (2a.. x.-b.) 
^  j=IlJ  J  1 


a  u. 


(64) 


In  matrix  form 


>  ww'x  -  ^cw  +  A' Ax  -  A»b  a  0 


(65) 
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where  w*  «  (w^, .  *  w^,)  is  the  vector  of  quadrature  weights.  Solving 
for  the  vector  x. 


x  r  (A1  A  +  y^w1 ) 


(A«b  +  ^cw). 


(66) 


If  we  define  the  unit  vector  e'  »  (1, ...  1),  then  the  vector  form  of 
the  area  constraint  solution  is 


x  a  (A* A  +  ^ee')”*  +  xce)  (67) 

where  w.al,  i  a  1,  •  •  • ,  M. 

Evaluating  P-,,(x)  by  substitution  of  (67)  into  (63)  we  obtain  the 
M 

following  quadratic  form  for  the  error: 

e  *  b*  {  I  +  A[  A* A  +  7vee»  ]"*  A*}  b 
+  2b{  t^A[  A’A  +  ^ee*  e}  c  (68) 

+(  ^  +  ^e'rA*A  +  ee*  j"1  e)  c^. 


This  may  be  written  in  the  form, 
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*  *  h'%-.h  +  2b'PMc  +  vc~  (69) 

which  is  referred  to  in  the  dynamic  programming  solution. 

To  begin  the  derivation  of  'he  dynamic  programming  solution,  we 
specify  that  1  <M<N.  For  M  =  1  there  results  the  simple  problem  of 
minimising 


N 


RjWntYfC)  +  s  (a^-h.)  . 


(70) 


The  minimizing  value  of  x^  is 


v 


bCwi  9  b'au; 

XWi2  +  ..(1)a(1) 


(71) 


where  a<(1)  =  (au>  a^ . aN1)  and  in  genera!  a'(M)  • 

th 

and  represents  the  M  column  of  the  A  matrix.  (x)  evaluated  at  the 
minimizing  value  of  x^  may  be  written  in  the  form 


rain  Rj(x)  *  b'Qjb+  2b*  p^c  +  r^c 


(72) 


where 


q  =. 


1  "  2  » 
*W1  +a’(l)a(l) 


pl  = 


-XWia 


*W1  +a.(l)a(l) 


2  2 
Wj 

?1  *  — J 

^W1  +a\l,a<l) 


The  recurrence  relationship 


*w(bic)  *  niiu[fw  .(b-a._-.x-_,  c-w_  x.  _)] 
M'  xMl  M-l  (M)  M*  M.  M'J 


is  stated  for  M>  2,  with 


fjtb.c)  =  rain  Rj(x). 


This  recurrence  relation  is  derived  in  the  Appendix  using  the  optimality- 
principle. 

Referring  to  the  definitions  previously  discussed,  the  parameters  b 

and  c  represent  the  state  variables,  the  functions  {  f^(b,  c)}  correspond 

to  tt  ^  criterion  functions,  and  the  optimal  policy  corresponds  to  the 

choice  of  x.  or  transformations  {b-a..,x.,  c-w.x.) . 

i  (i)  i*  i  i 

The  general  form  of  the  error  for  the  vector  solution  (69)  is  similar 
to  fj(b* c)  shown  in  (74)  and  (72).  It  is  a  simple  matter  to  prove 
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inductively  that  the  general  form  of  f^(b,  c)  is 


fM(b,c)=b.QMb  +  2b>pMc  +  rMc  . 


(75) 


To  solve  for  the  general  scalar  value  x,  „  the  recurrence  relation 

_________ 

(73)  is  used  with  (75) 


2  i 

b'Q  b  +  2b»p_  .c  +  r.  ,c  =  min(  (b-a._  _,x.  .)  C.  ,  (b-a,_  r.x  ) 

M  rM  M  x,,  (M)  M  M-l  .  (M)  M 

m 

+  2<b  -  a(M)V '  PM-i<C-WMXM>  ,76) 


+  rM-l(c-WMrM>  1  ‘ 


Expanding  the  expression  on  the  right. 


b*Q  b  +  2b’p  c  +  r,..c^  a  min  (b’Q. .  .b  +  Zb’p..  .c  +  r  ,c^ 

M  M  M  x  M-l  rM-l  M-l 

M 

+  X  mJ  rM-lWM  +  2wMa'(M)t’M-I  +  (77* 


*  2xMtCWMrM-l  +  +  K 


For  convenience  the  first  expression  on  the  right  in  brackets  is  defined 
as  G  and  the  second  expression  in  brackets  is  defined  as  D,  Differentiating 
the  expression  on  the  right  with  respect  to  x^  and  equating  the  result 
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to  zero, 


r,,  s  —  a 
M  D 


G_  a  CVm-1  *  ^'(M^M-l  *  WMbtpM-l*at(M)QM-lb  (?8) 
rM-lWM  +  2WMa,(M)PM-l+a*(M)QM-la(M) 


All  of  the  quantities  in  this  expression  are  known  except  CL  ,  ,,  p_ . 

M-l  M-l 

and  r  ,,  but  recurrence  relations  for  these  quantities  can  be  established 
by  evaluating  f^(b,  c).  Substituting  (78)  into  (77), 

b»Q  b  +  2b*p  c  4  r  c  *  a  b‘Q  tb  +  2b*p  c  +  r.  ,c 
M  M  M  M-l  M-l  M-l 


Expanding  G  and  equating  quadratic  coefficients  the  following  recurrence 
relations  are  evident: 


Q  =  Q 
M  M-l 


(WMPM-1  +“(M))(WMPM-1  +  W' 


rM-l  WM  + 


2,Vm  +  “m 


WMPM-1  rM-lWM  +  PM^ 

PM  *  PM-1 - 2~ - “7 - 

rM-lWKi  +  2WMPM  +  kM 


rM  =  rM-l 


(rM-l  WM  +  PM> 
rM-lWM  +  2WMPM  +  kM 
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where  for  computational  purposes  we  have  defined  auxiliary  quantities 


“(M)  =  °M-1  a(M)' 

PM  S  P'm-1  a(M)' 

and 

*  a'(M)  °M-la(M)* 


Recalling  tne  expressions  for  p^,  and  r^  in  (72)  we  see  that  if 
a  lt  P0  s  0,  and  =  x»  then  the  relations  just  developed  (80)  hold 
for  1<M<N. 

If  the  stability  control  parameter  \  =  0,  then  (78)  reduces  to 


a,(M)°M-l 
M  a'(M)QM-la(M) 


(81) 


Since  the  scalars  r^  =  0,  i  =  1, . . . ,  M,  and  the  vectors  p^  =  0,  isl,,,M 
M,  then  for  ^  =  0  the  recurrence  relations  (80)  reduce  to 


g(Mf*(M) 


(82) 


for  1<M<N,  with  Q  =  I. 

"  “  o 

The  general  relations  (78)  and  (80)  constitute  the  dynamic  programmir^ 
solution  for  the  restoration  problem  using  the  area  constraint  ao 
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previously  described.  These  relations  along  with  the  relations  developed 
by  Bellman,  et  al,  (1965)  for  the  prior  vector  constraint  were  used  in 
actually  simulating  the  restoration  process.  The  results  of  these 
simulations  appear  in  a  later  section. 

The  major  advantage  of  dynamic  programming  over  the  matrix  inverse 
solutions  is  that  no  matrix  inverses  are  needed.  As  a  further  evaluation, 
an  apparent  limitation  is  discussed  below. 

Pseudo-time  and  dynamic  programming 

Time  as  used  in  dynamic  programming  is  usually  used  as  an  index 
to  indicate  the  availability  of  information.  In  applying  dynamic 
programming  to  the  restoration  problem,  time  in  this  sense  is  an 
artificial  index  since  in  most  cases  all  of  the  information  needed  to  effect 
a  solution  is  available  at  the  same  time.  An  elaboration  of  these  state¬ 
ments  follows. 

In  business  problems,  as  well  as  in  many  other  problems  because  of 
the  sequence  of  time,  only  limited  information  may  be  available  to  operate 
upon  at  a  certain  stage  of  the  problem,  and  in  these  cases  one  can  do  no 
better  than  act  on  the  available  information.  For  these  cases  the  dynamic 
programming  solution  is  optimum  in  the  sense  of  using  all  of  the  available 
information.  However,  the  stepwise  minimization  as  noted  in  (70)  is 
accomplished  as  if  only  the  first  column  of  the  A  matrix  is  available, 
when  actually  the  entire  A  matrix  was  available.  There  is  some  cost  for 
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this  neglect.  If  this  cost  does  not  compensate  for  the  computational  error 
involved  in  matrix  inversion,  then,  neglecting  other  factors,  one  would 
not  use  the  dynamic  programming  approach.  This  tradeoff  has  not  been 
investigated  in  its  entirety,  but  the  evidence  from  the  results  available 
indicates  that  the  dynamic  programming  solutions  are  not  as  well  euitod 
io  the  optics  problem  as  the  matrix  inverse  solutions. 

Combination  of  dynamic  programming  and  matrix  inversion 

As  a  continuation  of  the  above  discussion,  one  logically  may  ask  if 
there  is  a  procedure  one  can  use  to  combine  the  properties  of  dynamic 
programming  (which  enable  one  to  treat  large  dimensionality)  and  yet  u  je 
more  A  matrix  information  per  object  estimate.  The  author  has  developed 
a  generalization  of  the  dynamic  programming  procedure  which  allows  one 
to  restore  the  object  vector  as  groups  of  subvectors.  The  solution  for 
each  subvector  involves  the  inverse  of  an  increasingly  greater  submatrix 
of  the  original  A  matrix  and  thus  utilizes  more  available  information  than 
the  scalar  by  scalar  solution  previously  considered.  It  appears  that  this 
generalization  may  be  applied  to  all  of  the  dynamic  programming 
constraints  used  thus  far. 

Because  of  the  computer  programming  complexity  involved  and  the 
fact  that  the  matrix  inverse  procedures  were  computationally  sound  for 
the  A  matrix  dimensionality  considered,  this  generalization  was  not 
investigated  further.  However,  future  efforts  should  not  overlook  this 


'""BirniimHIilttlKHItimllKIJl'lUltCltSOwiliio: . . 


generalization,  particularly  in  regard  to  the  greater  A  matrix 
dimensionality  inherent  in  two  •dimensional  object  restoration. 
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ANALYTIC  RESULTS  WITH  NOISE  FOR  SPECIFIC  CASES 

In  addition  to  the  theoretical  implications  that  this  section  conveys, 
this  section  is  presented  in  order  to  emphasize  the  trade-off  between 
diffraction  and  noise  in  regaining  object  information  which  has  been 
corrupted  by  these  effects*  Two  rather  separate  looks  at  this  trade-off 
have  been  investigated*  In  each  case  different  amounts  of  a  priori 
information  about  the  object  are  assumed* 

The  first  case  involves  the  problem  of  discriminating  between  two 
different  objects*  It  is  known  a  priori  that  there  is  a  choice  between  one 
of  two  objects  and  further  that  the  objects  are  specified  a  priori  as  a 
one-point  source  and  as  a  two-point  source  object.  In  the  second  case 
it  is  assumed  that  the  object  consists  of  two  point  sources,  but  because 
of  diffraction  effects  the  separation  between  the  two  point  sources  is  not 
known*  For  the  second  case  the  relationship  between  the  error  variance 
of  an  estimate  of  the  separation  and  diffraction  is  shown.  It  is  important 
that  the  a  priori  information  used  in  each  case  is  emphasized  so  that 
we  can  better  evaluate,  qualitatively  at  least,  the  coBt  involved  in 
regaining  various  amounts  of  information  from  the  observed  image. 

A  second  point  to  be  mentioned  is  that  tba  results  presented  in  this 
section  are  largely  theoretical,  while  those  that  follow  treat  the  mors 
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practical  aspects  involved  in  the  general  restoration  problem  when  there 
is  virtually  no  a  priori  object  information* 

Detection  Error  vs.  Diffraction 

The  problem  of  optical  discrimination  or  detection  has  been  considered 
by  several  authors  (Helstrom,  1964;  Harris,  1964b;  Rushforth  and  Harris, 
1966; .  n,  1966)*  In  this  treatment  several  assumptions  are  mace 
which  enable  the  analytical  solution  to  the  detection  problem  to  be  obtained. 
This  is  done  to  present  a  rather  complete  example  of  diffraction  and 
noise  effects  in  the  detection  problem. 

The  noise  model  and  detection  scheme 

To  begin  the  problem,  we  assume  that  the  image  intensity  distribution 
is  due  to  one  of  two  known  objects  plus  background  noise.  We  consider 
only  the  one -dimensional  version  of  Figure  1  and  that  the  objects  are  the 
point  source  6(a)  located  at  the  origin  and  the  sum  of  two  point  sources 
adta-rjj)  +  ^(a-t^)  located  at  and  respectively.  The  images  due 
to  these  objects  are  denoted  respectively  p(|)  aud  q(£).  To  enable 
analytic  results  to  be  obtained,  the  point  spread  function  is  taken  as 

h(|)  s  -U  e'e  /2<r  .  (S3) 
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Thus  using  (27)  the  images  are 


1  -e2/2*r2 

P(£)  s  1—T  e 

(84) 

•j2mr  ' 

and 

Id),  -d? 
jzvd 

(85) 

With  this  formulation  the  two  states  of  nature  are 

b(4)  .  p($)  +  n(l) 

(86) 

«,*  b<$>  ■  «Cft)  +  *<4) 

(G7) 

where  n(£)  is  the  additive  background  noise.  The  block  diagram  of  state 

of  nature  is  shown  in  Figure  2. 

It  is  further  assumed  that  the  noise  vector  n  (the  discrete  version 

of  n(£))  is  white  and  Gaussian  with  mean  n  =  0  and  covariance  matrix 
2 

K  =  <r  I,  The  more  practical  Poisson  noise  model  used  with  a  detector 
n  n 

is  discussed  in  the  Object  Estimation  section.  Since  n  is  Gaussian,  the 
conditional  density  functions  of  the  image  are  also  Gaussian  and  axe 


f(b/w,)-N(p,<r  2J) 
x  n 

f(b/«2)-N(q,«rn2I) 

(p*>) 

51 


where  N(n,  K  )  indicates  a  joint  Gaussian  density  with  mean  vector  (i  and 
covariance  matrix  K^,  and  b,  p  and  q  are,  respectively,  the  discrete 
versions  of  b(£),  p(§  j  and  q(£). 

I 

Using  the  techniques  of  statistical  decision  theory  (Middleton,  I960), 
we  can  at  this  point  specify  the  optimum  discrimination  procedure.  If 

! 

the  two  states  of  nature  are  equally  likely  a  priori,  and  if  the  costs 
associated  with  the  two  types  of  error  (i«e.,  choosing  u.  when  w,  is  tnv% 

I 

and  vice  versa)  are  equal,  then  it  can  be  shown  that  the  optimum  decision 
scheme  (in  the  sense  of  minimising  the  average  cost  or,  in  this  case, 

V 

the  error  probability)  is  the  following: 

choose  w.  if  &(b)  >  1 

(90) 

choose  c*2  if  £(b)  <  1 

where  j 

(b)»f(b/Wl)/f(b/«2)  (91) 

is  the  likelihood  ratio  of  the  observed  vector  b.  Since  the  natural 
logarithm  is  a  monotonically  increasing  function  of  its  argument,  an 

£ 

equivalent  statement  is: 


l 


£2 


choose  c>.  if  L(b)  >  0 

(92) 

choose  «2  if  L(b)  <  C 

where 

L(b)afnt(h).  (93} 


Noting  that  f(b/u^)  and  ffb/w^)  are  Gaussian,  substitution  of  (88)  ard 
(89)  into  (93)  results  in  the  decision  procedure: 

choose  if  (p-q)b  >  ijq'q-P'p) 

choose  (i>2  otherwise*  (94) 

Equation  (94)  is  just  the  discrete  form  of  the  correlation  detector  or 
matched  filter,  matched  in  this  case  to  the  difference  between  the 
diffracted  objects* 

Equation  (94)  specifies  the  data  processing  necessary  to  make  an 
optimum  decision*  However,  we  wish  to  investigate  more  thoroughly  the 
consequences  of  using  this  procedure  as  diffraction  and  noise  vary.  To 
accomplish  this  we  evaluate  the  probability  of  error  associated  wi*h  ti  e 
discrimination  procedure  just  described. 
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Evaluation  of  the  die  crimination  procedure 

There  are  two  types  of  error  which  a  discrimination  procedure  can 
lead  to:  Choosing  when  «2  is  the  true  state  of  nature  (i.  e. ,  calling 
the  two  point  source  a  one  poin?  source  in  our  case),  and  choosing  w2  when 
is  true*  If  decisionprocedure  (94)  is  used,  the  error  probability  for 
the  equally  likely  objects  is 


Pe  *  iP(G>P/u»2)  +  iPiG^/Wj) 


(95) 


where  G  «  (p-q)»b  is  the  test  statistic,  p  a  i(q*q-p'p)  is  the  freehold, 
and  f(G/wj)  is  the  conditional  density  of  G  given  that  is  the  true  state 
of  nature.  To  evaluate  the  error  probability  we  must  find  the  conditional 
densities  f(G/«2)  and  f(G/u>j)c 

The  Gaussian  noise  model  allows  these  conditional  densities  to  be 
found  by  noting  that  the  scalar  statistic  G  is  a  linear  combination  of 
independent  Gaussian  random  variables,  and  thus  G  itself  is  Gaussianly 
distributed*  To  specify  completely  the  densities  f(G/w^),  we  need  only 
find  the  conditional  means  E(G/<*>^)  and  the  conditional  variances 
var(G/uj),  These  conditional  means  and  variances  are,  in  this  case. 


easily  found.  They  are: 
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H-j  =  EvG/w^  =  (p-q)'p 

«a2  aE(G/«2)  r  (p*q)'q 

2  2  1 

ffj  a  var  (G/w^)  a  <r  (p-q)  (p-q) 

2  2  1 

o"2  »  var  (G/«2)  a  «rn  (p-q)  (p-q).  (96) 


Using  (96),  together  with  the  Gaussian  probability  density  function, 
the  probability  of  error  is 


'if. 


«  -<u1i2*2/2<r22 


e  vhT 


du  +  i  “Taiga-  1  1  —  ■  du* 


“CO 


l2mr  1 


(97) 


Further  algebraic  manipulation  yields  the  simpler  form 


r  12 

P  «  /  * - -  du 


•  / 

d/2 


J2ir 


(98) 


where 


d 

2 


(99) 
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Equations  (98)  and  (99)  axe  the  important  relations*  In  the  work 
that  follows  several  quantities  are  defined  which  will  enable  us  to  extract 
a  good  deal  of  information  from  these  relations  about  the  probability  of 
error  as  the  diffraction  and  noise  vary* 

In  order  to  simplify  (99)  recall  that  (04)  and  (85)  are  the  continuous 
versions  of  p  and  q*  Using  a  suitable  definition  as  to  the  transformation 
of  continuous  quantities  to  discrete  and  vice  versa,  it  can  be  shown  that* 
when  p  and  q  are  transformed  back  to  their  continuous  versions  (99) 
becomes 


d 

2 


(ICO) 


At  this  point  note  that  the  quantity  in  braces  may  be  defined  as  the 
"energy -to-noise"  ratio 


00 

EUR  =  d2  =  -ij  f  [p(5)  .  q<5)]2  d{  (101) 

<r  y 
n  oo 

where  the  term  "energy"  refers  to  the  integral  square  of  the  difference 
image*  (The  usual  definition  of  this  quantity  as  it  applies  to  electronic 
signals  is  the  signal-to-noise  ratio  (SNR)*  However*  in  this  paper  we 
reserve  SNR  for  the  more  literal  car'  re  signal  refers  to  the  true 
observable  image  intensity  or  counts  and  noise  refers  to  the  sample  mean 
absolute  variation  of  the  noise  intensity  or  counts* )  * 
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fr 

m 


The  analytical  formulation  of  the  detection  problem  is  found  by 
substituting  the  images  (84)  and  (85)  into  (100)*  Thus  d/2  becomes 

-fru-ni)2  2  ?  7  7, 

3  +  .  ST-.  4*"*  +  ) 


(102) 


Note  that  d/2  depends  upon  both  and  and  not  merely  on  their 
separation  (which  we  denote  y  ■  as  might  have  been  predicted. 

Thus  d/2  and  P  both  depend  not  only  on  ~te  noise  variance  and  diffraction 

© 

but  also  on  the  a  priori  choice  of  m  and  t^. 

It  is  not  immediately  clear  from  (102)  just  how  diffraction  plays  •» 

role  in  d/2  and  P  .  In  order  to  bee  this  role  more  cle-  ly  consider  the 
© 

case  when 


Now  the  importance  of  the  diffraction  ratio  R  as  defined  in  (28) 
becomes  more  apparent  since  R  is  in  this  case 


far 

r  s  — 

y 


where  k  is  an  arbitrary  constant  which  together  with  <r  specifies  the 


(103) 


(ICi) 


1 


"width"  of  the  point  spread  function  (i*  e* ,  kr  controls  the  aperture 
opening)  and  y  being  the  separation  is  defined  as  the  source  width* 
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The  basic  trade-off  between  diffraction,  noise,  and  P^  can  be  seen 
rather  easily  now*  V/hen  R  becomes  large  (y  0  for  fixed  o ),  then,  as 
one  would  suspect  d/2  approaches  sero  whish  drives  the  P  to  l/2--tbe 
worst  case*  Alternatively,  when  R  becomes  small  (y  oo  for  fixed  or) 
then 


d 

2 


(105) 


and  the  probability  of  error  reduces  to  a  quantity  which  depends  only  upon 

2  2 

the  noise  variance  <r  *  Not  until  <r  *♦  0  does  P  approach  sero  in  this 

n  n  e 

last  case.  (The  ratio  R  may  also  be  varied  by  fixing  y  and  varying  <r* 
Howe  <  £,  as  discussed  previously,  die  qualitative  results  are  similar, 
and  we  have  considered  throughout  the  paper  a  fixed  optical  system  which 
views  various  objects*) 

In  addition  to  the  above  discussion  on  limiting  cases,  we  wish  to 
consider  a  more  detailed  presentation  of  the  variation  between  the 
pertinent  quantities*  Figures  3-6  comprise  the  resulting  presentation* 

In  order  to  vary  die  a  priori  information,  i*e*,  and  r^i  three  cases 
were  considered:  Case  1,  rjj  ■  ca8e  2,  rjj  *  -1/4  and  case  3, 
fjj  ®  0  and  *  y.  The  variation  of  R  was  appropriately  accomplished 
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by  fixing  the  aperture  (<r  a  1/2)  and  setting  k  a  4;  thus  the  point  spread 
function  (83)  contains  95%  of  its  area  within  the  “aperture  width"  of 
4<r  ■  2. 0. 

The  variation  of  the  error  probability  with  diffraction  is  shown  in 

Figure  3.  The  limiting  cases  discussed  previously  are  evident:  (1)  As 

R  becomes  smaller  then  P  for  all  of  the  curves  approaches  a  quantity 

© 

which  depends  upon  the  noise*  note  that  this  quantity  changes  noticeably 
2 

as  cr  for  case  1  decreases  from  1. 0  to  •  001;  (2)  for  R  large  then  P 

a*  e 

approaches  0. 5.  An  interesting  feature  shown  in  this  figure  is  the 

trade  -off  between  a  priori  information*  diffraction  and  for  constant 
2 

noise  (r^  •  1. 0).  For  R  small  case  1  (the  two  symmetrical  point  sources 

compared  with  a  source  at  the  origin)  has  the  lowest  P  ,  but  for  R  large 

© 

case  1  exhibits  the  highest  P  • 

© 

This  result  may  contradict  000*8  intuition  at  first  glance*  but  as 
the  images  are  studied  for  the  3  cases  it  may  be  visually  seen  that  for 
R  small  it  is  easier  to  discriminate  between  the  images  for  case  1  than 
those  for  case  3  or  case  2.  On  the  other  hand,  for  R  large  enough  that 

l 

tixe  two  symmetrical  sources  of  case  1  are  located  "under"  the  point 
spread  function  (y  <  2, 0)  then  it  is  easier  to  discriminate  between  the 
images  of  case  2  or  case  3  than  those  of  case  1  because  one  of  the 
asymmetrical  sources  of  case  2  and  case  3  remains  outside  of  the  point 
spread  function  width  longer.  That  is,  R  has  to  become  larger  (for 

| 


robability  of  error  vs,  the  reciprocal  of  the  diffraction  ratio  R. 
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case  1  than  case  2  or  case  3)  for  both  of  the  sources  in  case  2  or  case  3 
to  be  located  within  the  point  spread  function  width. 

2 

To  further  illustrate  the  effects  of  the  noina  level  <r^  on  P  Figure  4 
2 

shows  P  vs.  log  (l/<r  )  for  fixed  large  value  of  diffraction  (R  »  10. 0)* 

c  n 

This  figure  clearly  illustrates  that  even  for  R  large,  P  **0  for  sufficiently 

© 

2 

small  e  ,  thus  confirming  the  notion  that  noise  nr  t  diffraction  is  tee 

true  limitation.  It  is  also  interesting  to  note  the  advantage  of  a  priori 

2  2 

knowL  dge  in  that  the  figure  indicates  a  savings  of  about  10  in  or^  (ENR) 

when  case  3  is  compared  to  case  1  for  P  s  0, 0. 

© 

Figure  5  Indicates  the  db  loss  incurred  by  increasing  diffraction 

while  maintaining  a  referenced  at  R  *  0. 0.  This  figure  shows  that 

the  diffra  ction  effect  acts  essentially  like  a  noise  amplifier  and  indicates 

that,  regardless  of  the  initial  noise  level,  the  effects  of  diffraction  alone 

can  cause  considerable  ENR  loss.  For  example,  if  the  discrimination 

procedure  produce:!  a  satisfactory  P^  for  R  »  .4,  then  to  maintain  this 

same  P  for  R  *  4.  0  would  require  a  21, 8  db  increase  in  ENR. 

© 

In  reiteration  of  the  concept  of  diffraction,  a  small  R  does  not  imply 
that  Pe~*0  since  as  K"*0  the  P#  still  depends  upon  the  uncertainty  caused 
by  noise.  In  this  respect  the  for  R  small  could  still  be  rather  high 
(perhaps  even  •  5).  Thus  the  db  loss  rhown  in  Figure  5  is  that  loss 
incurred  or  similarly  the  increase  necessary  to  maintain  a  referenced 
o'JLy  with  respect  to  diffraction. 
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We  tlto  desire  to  indicate  the  variation  of  the  db  loss  incurred  when 
2 

the  noiee  (<tq  )  as  well  as  diffraction  varies*  Figure  6  illustrates  this 

? 

variation*  In  this  figure  the  ENR  is  referenced  by  fixing  both  R  and  <r^ 

* 

such  that  P  a  *  05*  Two  points  axe  illustrated  ’,y  this  figure:  (1)  We 
o 

2 

can,  as  was  also  illustrated  in  Figure  3,  trade  tr  for  diffraction  and  still 

n 

'  maintain  die  same  P  (this  is  seen  as  we  observe  for  the  0  db  reference 

o 

2 

.joints  that  by  decreasing  <tq  from  1. 0  to  •  001  we  can  increase  R  by  a 

factor  of  1 0  while  still  maintaining  a  P  a  *  05)*  (2)  The  large  diffraction 

0 

effects  are  rather  severe  and  could  impose  a  severe  limitation  to  th; 

detection  process  unless  the  noise  can  be  greatly  reduced* 

The  results  of  ie  detection  problem  just  considered  afford  a  rather 

complete  look  at  die  basic  problems  imposed  by  diffraction  and  of  the 

2 

basic  relationship  that  noise  (o  )  can  be  traded  for  diffraction  (R>*  The 

n 

next  section  asks  a  related  question  with  regard  to  the  information  to  be 
regained  from  the  image  and  provides  further  insight  as  to  the  rather 
broad  applicability  of  the  relationship  just  mentioned* 

Separation  Error  Bound  vs.  Diffraction 

In  this  section  we  investigate  the  error  variance  associated  with 
estimating  the  separation  of  two  point  sources  when  no  a  priori  information 
is  given  as  to  the  location  of  these  two  sources.  The  problem  is  formu¬ 
lated  as  follows* 

E 

| 


■mrnX  *SE7« . " 
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We  consider  the  image  due  to  two  point  sources  located  at  and 
and  determine  the  covariance  matrix  of  the  errors  in  estimating  r\ ^  and  ti2* 
We  then  define  as  a  new  random  variable  the  separation  y  »  jT^-nJ  and 
find  the  error  statistics  of  this  new  random  variable  from  those  of  and 

V 

The  model  considered  is 


Mi)  *  +  Mi) 


(106) 


where  the  observed  signal  is  b(£)  and  n(|)  is  Gaussian  white  noise.  For 
these  conditions  Swerling  (1964)  has  shown  that  the  limiting  elements  for 
the  inverse  of  die  error  covariance  matrix  of  and  tj^  are 


Bij 


±  r: 
n  ,/ 


br\. 


br\. 


d$. 


(107) 


Again  we  consider  the  case  when  the  point  spread  function  is 
specified  by  (83).  Thus  when  x(or)  *  6{<*-iij)  +  dfc-r^) 


9(1# 


r-(i-t |j)2/2rJ  ^,)2/2<r2 

<  e  +  e 


(108) 


For  conveyance  let  the  covariance  matrix  of  tj  and  rj*  be  K.  Then 

X  M 

from  (107) 


(109) 


o 


and 

(k"‘)12  .  (K‘*)21  .  s - [arV).  (110) 

4N  nUc® 
o 

From  (109)  and  (110)  var  (^)  and  cov  (t)^  )  can  be  determined. 

Then,  as  referred  to  previously,  we  define  y  «  ti^ti  j  and  determine 
the  error  variance  of  this  parameter  (the  separation)  from 


var  (y)  n  var  (tjj)  +  var  fa2)  -  2  cov  (t^,^)  (111) 

The  resulting  limiting  error  variance  of  the  separation  is 


V  *  4 

•  /l  -e*^  /4<r  .  (112) 


This  expression  may  be  regarded  as  the  lower  bound  variance  for  all 
of  the  estimates  of  the  separation  which  are  minimum  variance  unbiased 


estimates 
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We  seek  to  investigate  the  variation  of  this  error  variance  with 
diffraction  and  noise*  As  in  the  detection  example  we  fix  c  ■  1/2  and 
consider  that  R  *  2.0/7*  The  limiting  cases  are  comparable  to  those 
resulting  in  the  detection  example*  For  R  large  (Y"*0)  then  the  error 
variance  approaches  infinity  and  the  separation  between  the  point  sources 
cannot  be  regained.  When  R  is  small  CY~*to)  then 

N 

V**  *s/lr  (113) 

and  here  as  in  the  detection  problem  the  error  depends  upon  the  noise 

2 

variance  (in  this  case  N  /2  *  <r  )•  Not  until  N  •*0  does  the  error 

on  o 

variance  approach  zero. 

To  further  illustrate  the  trade-off  between  the  error  variance  and 
diffraction  for  a  fixed  noise  level  we  consider  Figure  7  which  shows  the 
db  loss  vs.  1/R  (the  separation  in  this  case).  The  db  loss  is  referenced 
for  R  a*  1.64.  This  reference  was  chosen  because  the  general  curve  of 
error  variance  vs.  R  has  a  weak  minimum  at  R  **  1.64.  The  existence 
of  this  minimum  appears  to  violate  the  general  relationship  that  the 
error  is  proportional  to  R,  and  as  yet  no  physical  or  mathematical 
reason  is  apparent  which  explains  the  existence  of  the  minimum.  However, 
once  R  >  i.64  the  same  general  relationship  shown  previously  for  the 
discrimination  problem  exists  between  V  and  R.  For  example,  suppose 


IMIWIW  '  in  m.it«  1 1  u.m.  . 


db  loss  =10  log 


. . 
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that  we  were  satisfied  with  the  error  variance  achieved  for  R  ■  2.0; 
then  in  order  to  maintain  this  error  variance  for  R  b  10  we  must  provide 
for  an  increase  of  10,8  db«  Thus,  as  Figure  5  illustrates,  the  effect  of 
increasing  R  (decreasing  the  object  sise)  is  essentially  one  of  amplifying 
the  existing  noise  level,  and  in  order  to  maintain  a  constant  error 
variance  we  must  increase  the  SNR,  Also  noted  is  the  rather  severe 
limitation  imposed  by  excessive  diffraction.  This  same  comment,  as 
we  shall  demonstrate,  applies  to  the  numerical  solution  when  restoring 
incoherent  objects. 

The  next  section  presents  and  discusses  several  noise  models  which 
are  applicable  to  the  restoration  problem  and  develops  the  minimum  MSE 
estimate,  which  is  used  as  the  basic  solution  in  this  paper  to  restore 
objects. 
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OBJECT  ESTIMATION 

The  sequence  of  this  eection  is  as  follows.  First  we  discuss  ard 
present  the  noise  models  associated  with  the  restoration  problem.  The 
basic  sources  of  noise  consider  are  the  additive  background  noise  and 
the  multipUcitative  detector  noise. 

Next  we  consider  using  the  functional  form  of  the  multiplicitative 
density  function  in  estimating  the  object  vector  from  the  Baysean  approach. 
This  approach  offers  some  insight  as  to  how  a  priori  object  information 
might  be  encoded,  but  as  will  bv.  discussed  two  basic  drawbacks  to  this 
approach  are  that  the  additive  noise  is  considered  as  a  non-random 
variable  and  that  the  use  of  a  priori  object  information  in  the  inverse 
operator  is  not  possible  for  the  mathematical  model  considered. 

From  the  Baysean  approach  we  move  on  to  the  minimum  MSE 
approach.  Here  we  consider  the  additive  noise  as  a  random  variable  and 
ire  able  to  utilise  a  priori  information  in  the  inverse  operator  in  a 
strikingly  similar  manner  to  that  shown  by  Phillips  and  Twomey.  With 
a  reasonable  assumption  for  practical  problems,  it  is  easily  demonstrated 
that  the  MSE  estimate,  when  multiplicative  effects  are  considered,  differs 
only  by  a  multiplicative  constant  from  the  case  when  these  effects  are 
neglected.  Since  both  of  the  noise  models  can  be  effectively  represented 
in  the  MSE  estimate  and  since  a  priori  information  can  be  profitably 
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encoded,  the  basic  tool  for  object  restoration  used  in  this  pap  or  is 
considered  to  be  the  MSE  estimate. 

Although  the  MSE  estimate  effectively  accounts  for  noise  and  uses 
a  priori  information,  it  is  not  the  ultimate  estimate.  As  an  indication 
of  one  of  the  areas  where  the  MSE  estimate  could  be  improved  upon,  we 
conclude  this  section  by  discussing  the  minimum  distance  estimate  for 
the  object  x(a). 


Noise  Models  ,  _  V 

-  % 

Figure  8  shows  the  block  diagram  of  the  imaging  pi  ocess  and 
indicates  the  sources  of  noise. 

The  additive  model 

The  general  term  additive  noise,  as  previously  used,  has  referred 
to  the  observed  quantity,  either  intensity  or  mean  counts,  which  is  present 
when  the  object  radiation  is  absent.  The  radiation  which  causes  this 
intensity  or  mean  count  fluctuation  is  often  referred  to  as  ‘background* 
radiation  and  is  that  radiation  which  accounts  for  stray  light  fluctuations 
as  well  as  the  background  Intensity  in  which  the  object  intensity  distribution 
is  immersed.  We  assume  here  as  Helstrom  (1964)  does  that  this  noise 
is  basically  additive.  In  discrete  form  the  additive  model  is 


b  a  q  +  n  ■  Ax  +  n 


(114) 
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b  *  q  +  n  ■  Ax  +  n 


(114) 


where  n  is  the  additive  noise  vector. 

The  multiplicative  model 

When  the  radiation  strikes  the  detector,  which  is  often  a  photo¬ 
multiplier  tube,  a  single  photoemissive  surface  or  photographic  film, 
the  general  nature  of  the  reaction  must  be  described  in  statistical  terms. 
To  more  accurately  describe  this  model  the  surface  of  the  detector  is 
considered  as  being  divided  into  N  cells  each  with  incremental  «.rea  dA 
small  enough  such  that  the  intensity  impingi  ig  upon  any  one  cell  is 
constant,  throughout  that  cell.  In  practice,  it  is  convenient  to  characterise 
the  detector  surface  by  a  scalar  parameter  known  as  the  quantum 
efficiency,  which  we  shall  denote  as  q.  Under  suitable  practical 
conditions  the  mean  number  of  counts  observed  in  a  tLrr-.a  interval  of 
length  r  in  each  cell  may  be  taken  as 


-  B  Hi.  b 
i  hv  i 


(115) 


where  b^  is  defined  as  the  i**1  incoming  intensity  component  due  to  signal 
plus  noise  n.,  h  is  Plank's  constant  and  v  is  the  mean  frequency  of  the 
quasichromatic  radiation.  Furthermore,  the  counts  (number  of  resulting 
photoelectrons)  can  be  as sumed  for  suitable  conditions  which  have  been 


. . .  n . . 
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verified  in  practice  (Iviaadel,  1959;  Goodman,  1965)  to  obey  the  Poisson 
distribution 


ptyty 


k  -Db 
(Dbt)  e 

V 


(116) 


where 
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D  B 


(117) 
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It  also  appears  to  be  reasonable  to  asuume  that  the  counts  in  each 
cell  are  statistically  independent  from  those  in  another  cell  (Helstrom, 
1964),  Thus 


k.  -Bb. 
N  (Db.)*  e  1 


(118) 


where  k  is  the  count  ve  :tor  with  components  X, 

Equations  (115),  (116),  and  (118)  constitute  the  various  stages  of  the 
detector  noise  model* 


Temporal  variation 

As  a  further  discussion  of  the  noise  models  we  elaborate  on  the 
assumption  that  is  made  with  regard  to  temporal  variation.  In  the 
Imaging  Equation  section  it  was  stated  that  the  scalar  quantity  as 
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described  by  (2)  specifies  the  electromagnetic  field  radiating  from  &  point 
in  the  object  plane  as  a  function  of  distance  and  time.  In  carrying  on 
with  the  derivation  of  the  imaging  equation  we  assumed  that  no  information 
regarding  time  variation  would  be  considered.  To  make  this  assumption 
more  explicit  we  state  that  the  reason  for  this  assumption  is  that 
instruments  for  measuring  intensity  generally  cannot  follow  the  instan¬ 
taneous  fluctuations  at  optical  frequencies  (Helstrom,  1964).  Thus,  in  the 
strict  sense  the  observed  intensities  (b  and  n)  and  the  estimated  intensity 
x  are  time  averages.  This  does  not  imply  that  the  observed  quantities 
b  and  n  are  fixed  and  known.  That  is,  we  consider  that  a  priori  information 
is  available  about  these  observable  quantities,  and  in  the  general  case  we 
consider  that  a  priori  information  may  be  available  about  the  object  x. 

In  feet,  in  practice  the  time  average  will  only  be  taken  over  a  finite 
observational  interval,  so  that  from  a  practical  point  of  view  b,  n  and  x 
are  only  sample  means  and  not  true  means.  This  consideration  allows  us 
to  consider  a  more  flexible  model  of  add.  tive  noise  variation  in  that  the 
mean  noise  intensity  n  may  vary  from  observation  interval  to  Observation 
interval  as  well  as  spatially.  In  the  practical  sense  then,  n  is  still  to 

be  considered  as  a  random  vector  with  covariance  matrix  K  in  which 

n 

any  one  or  all  of  the  terms  in  are  nonzero  and  reflect  our  uncertainty 
about  the  true  noise  vector.  Further  discussion  concerning  time 
variation  is  presented  after  the  derivation  of  the  MSE  estimate  near  the 
end  of  this  section. 


i 
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Discussion  of  the  Baysean  Approach 

The  motivation  for  discussing  this  approach  is  due  to  the  general 
insight  which  this  method  of  estimation  theory  has  in  the  past  provided 
for  problems  of  this  general  type,  particularly  with  regard  to  the 
profitable  use  of  a  priori  information. 

The  term  Baysean  as  used  here  refers  to  the  general  estimation 
schemes  that  have  been  developed  using  the  density  functions  in  Bayes1 
theorem.  For  our  problem  Bayes'  theorem  is. 


f(x/k) 


Ptk/x)  f(x) 
p£(k/x)f(*)dx 


(119) 


“oo 

where  f(x/k)  is  the  a  posteriori  density  function,  P(k/x)  the  likelihood 
function,  and  f(x)  the  a  priori  density  function. 

Three  estimates,  the  maximum  likelihood  (ML),  maximum 
a  posteriori,  and  Bayes'  may  in  principle  be  obtained  from  (119)  once 
P(k/x)  and  f(x)  are  specified.  We  first  present  the  derivation  of  the  ML 
estimate.  Next  we  discuss  a  problem  which  arises  as  we  attempt  to 
determine  the  Bayes'  and  maximum  a  posteriori  object  estimates.  We 
will  indicate  how  this  problem,  which  is  one  of  mathematical  complexity, 
can  be  alleviated  at  the  cost  of  losing  the  use  of  a  priori  information  in 
tne  inverse  operator.  In  conclusion  we  will  present  the  form  of  the 
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resulting  Bayes*  and  maximum  a  posteriori  estimates  and  discuss  the 
insight  it  provides  as  to  the  weighting  o£  a  priori  information. 

The  ML  estimate  is  defined  as  that  set  of  parameters  x  (which  for 
this  estimation  procedure  are  regarded  as  fixed  but  unknown),  such  that 
the  likelihood  function  P(k/x)  is  maximised  (Mood  and  Greybill,  1963). 

In  obtaining  this  estimate  as  well  as  those  which  follow  in  this  discussion 
the  additive  noise  is  assumed  to  be  fixed  and  known. 

Finding  the  object  estimate  which  maximises  P(k/x)  in  (118)  is 
accomplished  as  follows.  Since  the  natural  logarithm  is  a  monotonic 
function  of  its  argument,  we  can  equivalently  maximise 


L(x)  «  t  n  P(k/x). 


(120) 


Upon  differentiating  L(x)  with  respect  ao  an  arbitrary  component  o2  x, 
say  x^,  and  setting  the  result  equal  to  sero  we  obtain 


.  t  .  .  N  f  k. 
3LM  .  s  _ !_ 

Ul  \S  "l 
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im 


3  0 


(121) 


M 

where  q.  a  2  a..x..  Noting  that  this  equation  holds  for  all  m  =  1,  • .  • ,  M, 
n  j=l  13  3 

then 


vA  s  0 


(122) 
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Pi 


where  v  ii  an  N  x  1  vector  whoee  component  is 


v.  ■ 


i  «li  +  ai 


-D. 


(123) 


Clearly,  if  A  is  square  (M»N)  and  nonsingular,  then  the  estimate  x 
can  be  found  from 


v>0, 


(124) 


Since  we  have  M=N  equations  of  the  form 


M  k. 

-  A  1 

25  ~n 

j.i  «  j  D 


-  n. 


(125) 


then  we  can  solve  for  the  ML  estimate,  which  is 


x  =  A  *(k /D-n). 


(126) 


The  ML  estimate  does  not  provide  any  new  information  which  would  over¬ 
come  the  basic  problem  of  an  unstable  A  matrix  whose  inverse  operates 
on  a  noisy  observed  vector.  It  does,  however,  indicate  that  the  ML 
procedure  results  in  the  straightforward  inverse  operator  which  operates 


i 
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on  a  modified  count  vector  (modified  by  dividing  out  the  multiplicative 
detector  effect  and  subtracting  the  fixed  and  known  background  intensity  n). 

As  was  demonstrated  by  Phillips  (1962)  the  addition  of  a  priori 
information  does  alleviate  the  unstable  A  matrix  problems;  thus  we  seek 
to  investigate  the  Bayes*  and  maximum  a  posteriori  estimates  which  use 
the  a  priori  density  f(x).  The  Bayes*  estimate  considered  here  is 
obtained  by  minimising  the  quadratic  loss  function,  while  the  maximum 
a  posteriori  estimate  is  designated  as  the  vector  x  which  maximises  the 
a  posteriori  density  f(x/k).  The  success  in  obtaining  these  estimates, 

41  88  mathematical  tractability  is  concerned,  depends  largely  upon 

the  form  of  f(x)s  Herein  lies  the  problem,  referred  to  previously.  The 
form  of  f(x)  must  be  such  that  one  can  mathematically  solve  for  either 
the  a  posteriori  mean  in  die  Bayes*  case  or  determine  a  unique  vector  x 
which  will  maximise  f(x/k).  Thus  far  we  have  been  unable  to  specify  a 
nontrivial  f(x)  to  enable  these  estimates  to  be  determined.  The  apparent 
difficulty  is  due  to  die  functional  relationship  between  the  object  x  and 
the  image  Ax. 


We  have,  however,  been  able  to  specify  an  a  priori  density  for  the 
image  vector  b  which  enables  the  determination  of  the  optimum  estimate 
for  the  image  b.  This  density  is  the  joint  gamma  density,  which  is 


N  u.  u  -b.c. 

_HC*  ^  e 

is  i  r(ut) 


(127) 
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where  the  parameter*  c  and  u.  specify  the  mean  b.  =  —  and  the 

*  i  i  c. 

2  * 

variance  var  (b^)  «  u^/c^  •  Here  we  have  assumed  that  the  intensity 
is  independent  from  b^  for  all  i  and  j. 

Utilisation  of  f(b)  as  specified  above  with  the  Poisson  distribution 
allows  the  determination  of  the  a  posteriori  density  f(b/k),  which  is  again 
a  joint  gamma  distribution  with  parameters  D  +  and  k.  +  Uj  which 
replace,  respectively,  the  parameters  Cj  and  u.  in  (12?)*  With  the 
specification  of  f(b/k)  in  the  form  of  the  joint  gamma  distribution,  we 
can  solve  for  either  the  Bayes'  or  the  maximum  a  posteriori  estimate 

.A 

for  the  image  b.  Using  this  estimate  we  consider  that  a  reasonable, 

A 

although  suboptimum,  procedure  is  to  operate  on  b  to  obtain  an  estimate 
for  x. 

This  procedure  is  easily  accomplished  for  the  Bayes'  case.  Since 
the  Bayes'  estimate  is  the  a  posteriori  mean  we  have 

-  ki + "i 

bi-EtbAl  *  T5TT  •  (128> 

i 


Using  the  proecdure  just  dib.  ussed,  we  have  N  equations  of  the  form 


„  M 

b.  *  S  a..x.  +  n. 
1  >1  « 1  1 


Vui 

D  +  Cj* 


(7  29) 
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For  A  non  singular  and  (M*N)  then  the  suboptimum  estimate  for  x  is 

xmAml[k  +b  -n]  (130) 

c  c 

where  k  and  b  are,  respectively,  N  x  1  dimensional  vectors  with 
c  c 

components 


and 


k 


D  +  c. 


ci  bi 
D  +  c. 


with  b.  *  u./Cj,  The  suboptimum  maximum  a  posteriori  estimate  for  x 
is  found  in  a  similar  manner  as  that  demonstrated  in  obtaining  the  ML 
estimate.  Its  form  is  essentially  that  of  the  Bayes*  estimate  shown  in 
(130). 

In  order  to  discuss  the  insight  Bayes*  estimate  provides  in  using  a 
priori  information,  we  consider 


-  N  -1 
x  *  S  a 
1  jel 


-S-fiV  £ . 

D  +  c,\D  /  D  +  c. 


j  “j 


(131) 


where  the  a^’s  are  the  elements  of  A*1. 
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Recall  that  the  constant  D  is  defined  in  (117)  and  contains  the 

observation  period  t.  If  t  become*  large,  D  becomes  large  and  the 

observed  counts  k.  approach  the  moan  counts  a.  and  thus  k  /D  approaches 
j  j  j 

the  mean  intensity  b  , 

<1 

The  mean  and  variance  of  the  j  ^  prior  gamma  distribution  are 

i  2 

Uj  j  J  Cj  *  ^u#>  generally  speaking  small  values  correspond 

to  large  prior  uncertainty,  since  the  variance  is  large  and  the  a  priori 

distribution  is  rather  broad  and  flat. 

Now  refer  to  (131)  and  consider  the  case  when  c .  is  small  and  D  is 

J 

large.  Here  we  see  that  the  estimation  procedure  weights  the  observation 

much  more  than  the  a  priori  mean  component  £.,  This  case  corresponds 

J 

to  very  large  a  priori  uncertainty,  and  naturally  we  have  little  confidence 
in  the  a  priori  mean. 

Alternatively,  consider  the  case  when  is  large  and  D  is  small. 

The  situation  is  the  reverse  of  that  just  considered.  We  weight  the  a 

«e 

priori  mean  b  heavily  now  and  almost  disregard  the  observation  k., 

1  J 

Here  the  estimation  procedure  reduces  to  one  of  essentially  "cleaning  up" 
or  enhancing  the  a  priori  mean 


Minimum  MSE  Approach 


In  this  section  the  estimation  problem  is  approached  from  a 
different  point  of  view,  which  does  not  require  the  specific  functional 


. . .  •'tlliaillllllimilllllllilillllllllllll . II . . 
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form  of  the  a  priori  density  or  the  likelihood  function  to  be  specified.  We 

now  assume  that  the  additive  noise  vector  ;  §  random  and  treat  two  basic 

types  of  MSE  estimation.  The  first  is  the  "classical"  least  squares 

approach,  which  assumes  that  no  a  priori  information  about  the  object  is 

available  and  that  the  additive  noise  has  aero  mean  and  covariance  matrix 
2 

°n  **  Secondly,  we  consider  the  more  general  case  when  the  noise 

covariance  matrix  is  K^,  and  we  take  advantage  of  encoding  a  priori 

object  information  in  the  form  of  the  prior  mean  x  and  the  covariance 

matrix  K  , 
x 

Before  beginning  the  two  approaches  the  noise  model  is  more 
explicitly  seated.  Recall  (115)  and  note  that 

M 

^/D.Sa  x  +nt  (132) 

jsl  1  J 
or 

bi=SaijXJ+ni*  (133) 

Since  n^  is  a  random  variable,  the  quantity  z^/D  or  b^  is  now  a  random 
variable.  If  n^  has  mean  zero,  which  we  cam  assume  without  loss  of 
generality,  then  the  model  is 


J . MEM . .  .'jlHMIWIViB 
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b  «b  +  n 


(134) 


where  b  a  Ax  and  £  is  the  mean  of  b. 

Classical  least  squares 

The  classical  least  squares  approach  presented  here  has  been  treated 
by  several  authors.  The  reader  is  refe*  red  to  Mood  and  Greybill  (1963) 
and  Deutsch  (1965)  for  turther  discussion*  It  is  instructive  to  mention 
that  in  the  general  case  treated  by  Deutsch  (1965),  the  assumption  of  a 
linear  relationship  may  be  nece^-ary  between  the  unknown  parameter  x 
and  the  obse  *ved  random  vector  b*  However,  in  the  restoration  problem, 
this  assumption  is  not  necessary  since  the  imaging  process  is  inherently 
linear. 

2 

We  consider  here  the  case  for  K  sir  L  The  parameter  x  is  chosen 

n  n 

to  minimize 


R(x)  =  (Ax-b)  (Ax-b) 


(135) 


or  alternatively 


N  M 


R(x)  n  S  (S  a  x -b.)  . 
islial  J  1 


(136) 
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The  parameter  x  is  determined  in  the  usual  manner  by  differentiating 

R(x)  with  respect  to  an  arbitrary  element  of  x,  say  x  ,  and  setting  the 

m 

result  equal  to  aero.  In  summation  form 


N  M  A 

S  a  (S  a  x  -b  )  *  0,  (137) 

i=l  1111  j=l  «  3  1 


since  this  expression  must  hold  for  all  m  *  1, . , , ,  M, 

A«Ax  =  A’b.  (138) 


If  (A'A)  is  nonsingular,  then  the  estimate  is 


vc  «  (AfA)"lA*b# 


(139) 


When  A  is  square  and  nonsingular 


a  .  -1. 

X  a  A  b. 


(140) 


Thus,  the  classical  least  squares  approach  results  in  the  same  basic 
inversion  form  as  shown  in  the  ML  case  and  as  seated  in  the  Invvging 
Equation  section*  Even  though  statistical  methods  have  been  used  in 
attaining  this  basic  form,  we  still  have  not  been  able  to  gain  much 
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additional  insight  which  would  enable  us  to  overcome  the  problems  of 
instability  and  noise  referred  to  earlier*  The  inclusion  of  a  priori 
information  as  shown  next  does,  however,  provide  considerable  insight 
into  these  problems* 

The  general  minimum  MSE  estimate 

Now  we  assume  that  a  priori  information  is  available  in  the  form  of 

the  prior  object  mean  x  and  the  covariance  matrix  K^,  Information  about 

the  noise  ve'  tor  n  is  assumed  to  be  available  as  -needed  in  the  covariance 

matrix  K  ,  and  we  assume  that  the  noise  mean  is  zero.  We  further  assume 
n* 

that  x  and  n  are  independent*  The  procedure  which  follows  in  obtaining 
the  general  result  (149)  was  first  applied  to  the  restoration  problem  by 
Rushforth  (1965)* 

We  postulate  a  linear  model  which  enables  the  use  of  K^,  x  and 
We  assume  that  x  and  the  random  intensity  vector  b  are  related  by 


x  =  Hb  +  g 


(141) 


where  H  is  an  MxN  matrix  and  g  is  an  Mxl  column  vector.  We  seek  H 
and  g  such  that  the  quantity 


MSE  =  E  [(x-x)«  (x-x)J 


(142) 
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is  minimized  where  x  represents  the  true  object.  Once  H  and  g  are 
determined,  then  the  optimum  estimate  x  is  found  by  substituting  these 
quantities  into  the  expression  for  x  above. 

Recall  (114)  that  b  a  Ax  +  n.  It  can  be  shown  (Rushforth,  1985)  that 
the  vector  g  can  be  written 


g  »  (I-HA)x. 


(143) 


Substitution  of  (141)  and  (143)  into  (142)  yields 

MSE  a  tr  H(Ax  +  n)  +  (I-HA)x-xJ 

[  H(Ax  +  n)  +  (I-HA)x-x]  '  j-  (144) 

where  tr  denotes  the  trace  of  a  matrix.  Upon  expanding  (144)  and 
performing  the  expectation  over  x  and  n  we  can  write  (144)  in  the  following 
form: 

MSE  a  tr/k[A  KA'  +  K  JH*-K  A»H*-HAK  +  K  (145) 

L  1  x  n4x  xxj  ' 

To  solve  for  H  we  complete  the  square  in  (145)  and  obtain 


88 


e  a  tr  /[ H(A K  A*  +  K  )*  -  K  A«(AK  A«  +  K  )"*] 

(  x  n  x  '  x  n  J 

[H(AK  A*  +  K  )^  -  K  A«(AK  A*  +  K  )“^]« 

*■  '  x  n7  x  x  n  J 

+  Kx-KxAl(A  KA>  +  Kn)‘*  A  Kj.  (146) 

Now  the  minimising  value  of  H  is  easily  determined  by  equating  either 
one  of  the  two  expressions  in  brackets  above  equal  to  zero.  Thus, 

H  a  K  A*(AK  A*  +  K  f1.  (147) 

x  '  x  n7 


Further  manipulation  reveals  that  H  may  be  represented  in  the  alternate 
form 

H  a  (A’K  "*A  +  K  “1)“1  A*K  *l.  (148) 

n  x  n 

Using  equations  (148)  and  (143)  for  H  and  g  the  general  expression  for 
x  is 

x  a  (A«K  _1A  +  K  “VW  -1b  +  K  _1x).  (149 

'  n  x  n  x 

An  important  quantity  to  consider  in  the  evaluation  of  how  well  x 


performs  is  the  MSE  which  results  when  x  is  used.  The  resulting  MSE 


89 


is  found  as  follows*  By  substituting  H  into  (146)  and  noting  the  equality 
of  (147)  and  (148),  the  MSE  expression  becomes 


MSE  e  trTlC  -(A»K  “*A  +  K  ml)~lA*K  -1 
L  x  n  x  ’  n 


AK 


1 

=]• 


(150) 


Thus, 

MSE  a  tr  J(A*K  “lA  +  K  "*1)"*1[(A,K  "*A  +  K-1)K  -A’K”1*  K  ll.  (151) 

(  n  x  n  x  x  a  x ) 

The  quantity  in  brackets  in  (151)  reduces  to  I;  thus,  the  general  expression 
for  the  MSE  reduces  to 


MSE  &  tr 


f(A'K, 


-‘a+k-v1]. 

n  x  '  j 


(152) 


Multiplicative  effects  for  the  high  SNR  case 

We  assume  in  this  paper  that  the  observed  quantity  is  either  intensity 
or  mean  counts.  That  is,  the  observation  interval  is  assumed  to  be 
long  enough  so  that  when  the  detector  is  present  the  actual  number  of 
counts  k.  is  a  very  good  approximation  for  the  true  mean  counts  z..  Thu? 
when 


k.« 

i 


zi* 


(153) 
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the  optimum  estimate  for  x  is 

x  .  [A>K‘lA  +  Kx‘1]*1tA'K^1(*/D)  +  (154) 

In  practice  when  deviates  significantly  from  z^,  we  must  consider 
a  more  sophisticated  model  such  as  that  treated  by  Austin  (1966)  in  which 
he  assumed  that  we  observe  counts  as  they  emerge  from  the  detector  but 
estimate  the  mean  intensity  x  of  the  object*  This  treatment  adds  insight 
(in  the  limiting  case  for  T  **  co  the  estimates  are  the  same),  but  indicates 
that  much  more  complexity  is  involved  in  the  restoration  process  when 
counts  are  observed* 

An  increase  in  T  essentially  implies  an  increr.se  in  SNR.  The 
evidence  presented  in  the  Analytic  Results  with  Noise  for  Specific  Cases 
section  and  the  Simulated  Object  Restoration  section  indicates  that  for 
reasonable  restoration  (for  R  >  1*0)  we  must  consider  the  high  SNR  case. 
These  results  indicate  that  the  SNR  must  be  high  enough  so  that  we  can 
neglect  the  more  complex  model  and  just  consider  the  general  case  (154;, 
By  high  SNR  we  are  implying,  generally  speaking,  that  for  the  diffraction 
range  considered  (R  >1*0)  the  SNR  must  be  high  enough  to  neglect  the 
detector  except  for  a  multiplicative  constant. 

Another  way  of  stating  the  high  SNR  assumptior  is  to  Mate  that  the 
SNR  is  high  enough  so  that  the  central  limit  theorem  applies  (Parzen,  I960) 


.  . . . . . . mm . nil.  . . .  wi*nn» . .  j.d;.uiMi>s.«wii.diaHuiau>daaiHi«>idlMiiiiijniMMri«MMiMiliilii|i*jiliaHl4H<«il>h-99i^KiilniMji!:-mliil:is:niAr<iiuliiiiiilnii]iuiiiiiailii« . . . MM . . 
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and  we  have  Gaussian  rather  than  Poisson  statistics.  In  this  case  it  is 

easily  demonstrated  that  if  then  the  use  of  Gaussian  distribution 

functions  for  f(z/x)  and  f(x)  result  in  a  Bayes*  estimate  or  a  maximum 

a  posteriori  estimate  that  is  equivalent  to  (154).  Thus,  the  assumption 

of  high  SNR  ties  the  Baysian  and  MSE  approaches  together. 

To  follow  up  the  earlier  discussion  on  temporal  variation,  the  high 

SNR  assumption  implies  that  the  observation  interval  is  long  enough  so 

that  n,  b  and  x  are  sample  means  with  respect  to  time  that  are  near,  but 

not  equv  1,  to  the  true  temporal  mean.  Thus,  K  and  K  reflect  our 

X  & 

ignorance  of  the  true  quantities.  As  is  the  usual  case,  large  elements 
in  K  and  K  imply  a  greater  uncertainty  than  small  elements, 

U  A 

Next  we  consider  special  cases  of  the  general  estimate  (149)  and  the 
MSE  (152)  and  din  cuss  the  trade-off  between  a  priori  information  and 
noise. 

Discussion  of  the  MSE  estimates— special  cases 

The  trade-off  between  a  priori  information  and  noise  is  easily  seen 
2  2 

when  K  s«r  I  and  K  »<r  L  The  general  estimate  (149)  becomes 
xx  &  n 


(15J; 


x  =  (A*A  +0"  ^1)  *(A'b  +  <r  ^/<r  ^ x), 

n  x  '  '  n  x  1 
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Consider  the  case  when<rx  >><rn  *  This  case  corresponds  to  large 
&  priori  uncertainty  ’T*re,  as  before  in  the  discussion  of  the  Bayes  and 


maximum  a  posteriori  estimates!  we  tend  to  weight  the  observation  vector 

b  much  more  heavily  than  the  prior  mean  x.  On  the  other  hand,  when 

2  2  ml 

<rn  >><TX  »  w®  a  noisy  observation  vector  b  and  we  weight  x  heavily. 

Expressions  (149)  and  (155)  gives  us  more  insight  than  (130)  in  that 


the  basic  inversion  operator  is  also  modified  in  accordance  with  a  priori 


information.  Note  the  similarity  in  form  of  (155)  and  solution  (47).  The 


analogy  between  the  various  quantities  in  these  two  equations  is  also 


evident.  Particularly  is  this  so  when  the  covariance  matrix  is 
compared  with  H  and  the  a  priori  mean  r  is  compared  with  the  a  priori 

! 

§ 

vector  p.  The  fact  that  Doth  of  these  approaches  result  in  the  same  basic 

2 

form  further  enhances  the  use  of  a  priori  information  in  the  inverse 

I 

i 

operator,  and  it  will  be  demonstrated  that  an  essential  part  of  the 


restoration  process  is  the  stability  control  offered  by  use  of  this 


information. 


In  order  to  show  how  the  general  case  (149)  reduces  to  the  classical 


MSE  estimate  (139)  and  to  introduce  the  problem  of  measurement 


selection,  consider  the  case  when  there  is  essentially  no  contribution  due 


to  a  priori  information.  In  this  case  the  estimate  (149)  may  be  written 


a  .1  .1  -1 

x  h  (A«K  A)  (A»K  b). 
n  n 


. . . 
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This  is  known  as  the  Markov  estimate  and  corresponds  to  the  minimization 
of  the  weighted  least  squares 


R(x)  *  (Ax-b)  K  (Ax-b) 
n 


(15?) 


where  as  before  Kq  is  the  noise  covariance  matrix  (Deutch,  1965).  When 
2 

Kn  a  <rQ  1,  the  white  noise  case,  the  estimate  in  expression  (156)  reduces 

as  it  should  to  the  classical  least  squares  estimate  of  (139)* 

•*1  »1 
When  the  elements  in  K  are  small  compared  with  those  in  A’K  A 

x  n 

the  error  (152)  reduces  to 


MSE  atrftA'K^A)"1], 


(158) 


When  Kao-  I  we  have 
n  n 


MSE  *«rn2  tr[(A'A)“l], 


(1 


KOI 

•»  /  } 


or  tiie  equivalent  form 


2  N  |  2 

MSE  » <r  E  l/x7 
n  ui  »  rti 


(16C, 


i  ax*-,.  ‘fess s2 


.:(«***  ««i<>  iwr-ar-  — irlHPt-milijl 
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where  the  are  the  eigenvalue*  of  the  A  matrix.  This  fcrm  (160)  was 
also  suggested  by  Barnes  (1966b)  as  an  important  form  to  consider  when 
studying  the  effects  of  diffraction  when  using  the  classical  inversion 
proce  s. 

When  K&  is  specified  in  (1 58),  only  the  A  matrix  can  be  varied  to 
reduce  the  error.  This  is  both  a  nontrivial  and  important  problem.  It 
is  nontrivial  because  it  is  very  difficult  to  solve  for  optimal  conditions 
that  will  specify  the  A  matrix  for  each  restoration.  It  is  important 
because,  as  will  be  shown,  the  MSE  can  be  reduced  considerably  (by 
several  orders  of  magnitude)  by  just  knowing  where  to  measure  the  image 
and  where  to  predict  the  object. 

*  Austin  (1966)  has  been  able  to  solve  for  the  optimal  condition  in  a 

special  case.  The  special  case  is  one  in  which  the  image  is  assumed  to 

be  the  summation  of  equally  spaced  diffraction  patterns  from  known  point 

2 

sources  which  are  separated  by  the  Rayleigh  distance  and  I.  For 

this  special  case  he  showed  that  the  optimum  image  measurement 
locations  are  above  the  known  point  source  locations  which  are  transferred 
to  the  image  plane.  In  this  case  A:I,  and  the  MSE  becomes 

MSE=<r2M.  (161 ) 

n 


i 


;! 
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In  this  paper  we  consider  the  more  general  and  more  important  case 
since  we  consider  here  that  the  object  is  to  be  both  unknown  and  continuous 
and  that  the  diffraction  effects  are  severe  enough  so  th  .t  we  must  predict 
at  points  separated  by  a  distance  which,  for  most  cases,  is  only  a  small 
fraction  of  the  Rayleigh  distance*  To  further  specify  the  case  we  consider 
in  this  paper,  we  reiterate  that  the  value  of  R  is  to  be  greater  than  unity. 
(The  case  Austin  (1966)  considered  holds  for  R  <1.0*  For  R  *  1*  0  the 
number  of  point  sources  M  *  3*  0, )  ForR  >1.0,  the  optimum  scheme 
for  choosin;  the  A  matrix  parameters  has  not  bee  solved.  However,  in 
order  to  provide  a  rather  complete  computational  procedure,  the  variation 
of  the  parameters  which  determine  the  A  matrix,  and  thus  the  MSE  hac 
been  extensively  studied  by  numerical  methods.  Results  of  this  study 
and  the  computational  insight  which  they  provide  for  the  working 
restoration  procedure  are  presented  in  the  secuon  entitled  “MSE 
Variation- -Choosing  the  A  Matrix  Parameters. " 

Perhaps  the  outstanding  problem  which  remains  to  bo  solved  is  that 
of  initially  estimating  the  size  of  the  diffracted  object.  This  important 
problem  is  discussed  further  below. 

Minimum  Distance  Estimation 

Minimum  MSE  estimation  considers  the  minimization  of  functional 
differences  between  the  predicted  object  and  the  true  object.  The  integral 
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equation,  however,  requires  that  we  also  choose  other  parameters  before 
the  vector  x  can  be  estimated*  For  example,  M,  N,  the  cri's,  the  £i's 
and  the  must  r U  be  chosen  before  estimating  ^*  Evidence  for  making 
intelligent  choices  for  all  of  these  parameters  is  discussed  later*  Ho  re 
we  present  a  minimum  distance  criterion  which  accounts  for  both  the  size 
and  functional  shape  of  the  object* 

Visually,  diffracted  objects  appear  smeared  or  spread  out.  The 

a 

size  of  the  true  object  is  not  known*  To  effect  a  better  restoration  than 
the  MSL  estimate  we  consider  choosing  a  and  x  such  that  the  following 
mean  square  distance  is  a  minimum: 


d2  a  E  ^[x(a  )-x(o )] 1  [x(a )-x(<*  )] 

+[£-c]i  (162) 


The  minimization  of  this  expression  implies  choosing  a  and  x  together. 

No  analytic  results  have  been  obtained  for  accomplishing  this  task*  The 

l 

problem  thusfar  has  been  that  the  size  and  object  functional  values  are  so 
intimately  related  that  we  cannot  estimate  one  without  knowledge  of  the 
other. 

In  the  absence  of  an  optimum  size  estimate,  the  object  size  and 
likewise  the  diffraction  have  been  initially  estimated  by  subtracting  the 
poiuw  spread  size  from  the  imag'd  size*  Using  this  initial  aize  estimate, 
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we  have  then  sequentially  estimated  the  object  size  and  shape  together* 
However,  it  appears  that  we  can  do  better,  and  future  efforts  undoubtedly 
should  consider  this  problem* 
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SIMULATED  OBJECT  REST  ORATIONS 

In  this  section  we  present  the  results  obtained  from  actual  numerical 
computations  performed  on  a  digital  computer  which  simulate  the 
various  restoration  solutions  developed  previously*  The  major  computer 
program  used  in  performing  these  simulations  is  presented  in  the 
Appendix*  All  of  the  solutions  which  utilise  a  priori  information  in  the 
inversion  process  have  been  investigated*  The  MSE  solution  with  the 
additional  computational  feature  of  iteration  has  been  used  most 
extensively* 

For  all  of  the  simulations  we  have  considered  the  one  -dimensional 
incoherent  point  spread  function 

hU-a)  *  sinc2(£-or)  (163) 


where  we  have  fixed  the  aperture,  as  previously  discussed,  so  that  the 
Rayleigh  separation  is  unity*  (Equation  (163)  describes  an  optical 
configuration  with  unity  magnification.  This  involves  no  loss  in  generality 
since  it  is  always  possible  to  normalize  the  object-image  coordinates 
so  that  the  magnification  is  unity.)  Using  (163)  the  imaging  equation  is 


m -f 


sine  (£-a)x(a)d a, 


-a 


(164) 
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Results  are  presented  on  the  interrelationship  between  the  diffraction 
level  (R),  the  noise  level  in  the  image  and  point  spread  function,  the 
choice  of  parameters  which  govern  the  A  matrix,  and  the  amount  and  type 
of  a  priori  information  used.  Also  considered  is  a  method  of  smoothing 
the  image  when  we  encounter  excessive  noise. 

In  this  section  the  results  of  the  interrelationship  just  mentioned  are 
discussed  as  the  section  proceeds.  For  a  summary  of  these  inter¬ 
relationships  and  the  computational  scheme  which  has  been  developed 
the  reader  is  referred  to  the  Summary  and  Conclusions  section.  The 
computational  scheme  is  also  discussed  in  the  next  section  entitled  "MSE 
Variation— ‘Choosing  the  A  Matrix  Parameters. 11 

Before  presenting  the  simulations,  it  is  necessary  to  present  a 
clearer  picture  of  the  simulated  restoration  process  by  discussing  three 
general  areas.  These  areas  consider  the  objects  to  be  restored,  the 
image  and  A  matrix  accuracy,  and  the  choice  of  a  restoration  improvement 
criterion. 


Objects  to  be  Restored 

Basically  four  objects  have  been  considered.  These  are  the  uniform 


pulse. 


,  .  /l.  0  for  -a<ar'<  a) 

*<“>  "  |  0  elsewhere*  ) 


(165) 
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the  smooth  pulse, 


x(ct) 


t 


cos  kiror  for  -a  <a  <  a 


0  elsewhere 


-1  _  a  ( 


J 


(166) 


two  uniform  pulses, 


x(or)  = 


*1,0  for  -a^  <a  <  aj 


V.. 


“a2  la2 

0  elsewhere 


(167) 


and  two  smooth  pulses. 


x(o) 


fcoa  krr(a-<p)  for  -a^  ^ar  <al  ) 
=  <!  *  -a?  <«  <a7  { 

-  2J 


0  elsewhere. 


(168) 


The  motivation  for  considering  these  basic  objects  was  as  follows. 

First,  recall  that  the  solutions  investigated  utilize  a  smoothing  parameter 

2 

which  in  the  MSE  case  depends  upon  the  noise  variance  <r ^  and  the  a  priori 
variance  <r  Phillips  <1962)  has  mentioned  that  his  solution  should  work 
for  x{a)'8  that  are  smooth.  In  order  to  gain  an  estimate  of  the  loss 
involved  when  the  objects  are  not  smooth,  various  smooth  (cos  knar)  and 
non-smooth  (the  uniform  pulse)  objects  have  been  chosen.  Secondly,  it 
is  interesting  to  determine  how  well  the  restoration  procedure  can 
restore  split  sources.  Such  objects  are  certainly  difficult  to  restore 
(perhaps  the  "worst  case"  would  be  a  large  number  of  split  sources 
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which  are  very  closely  spaced),  and  they  provide  a  more  absolute 
criterion  as  to  how  well  the  restoration  procedure  performs.  Furthermore, 
the  resolution  of  split  sources  which  have  been  excessively  diffracted 
provides  a  good  comparison  with  the  classical  Rayleigh  criterion. 

Image  and  A  Matrix  Accuracy 

As  has  been  previously  emphasised,  noise  is  the  limiting  factor  in 
the  restoration  of  optical  objects  of  finite  extent.  Thus,  in  the  simulations 
which  are  to  be  presented  it  is  important  that  we  specify  the  procedures 
used  to  obtain  the  image  and  A  matrix  and  to  provide  measures  as  to  the 
accuracy  of  these  quantities.  This  section  presents  this  information. 

Obtaining  the  images 

The  images  used  for  the  restoration  process  were  obtained  by 
numerical  integration  using  Simpson's  quadrature.  The  motivation  for 
using  numerical  methods  was  essentially  twofold.  The  integrations  for 
a  variety  of  sources  one  may  wish  to  consider  cannot  always  be  performed 
analytically.  Also,  since  the  restoration  process  is  one  of  "inverse- 
integration,  "  then  it  is  natural  that  one  would  learn  from  the  forward 
process,  particularly  in  regard  to  the  choice  of  the  number  of  points  for 
a  desired  accuracy.  A  third  consideration  is  that  a  numerical  procedure 
allows  one  to  determine  the  image  for  arbitrary  arguments,  while  the 
analytic  image  (if  available)  may  be  difficult  to  obtain  for  some  arguments. 


. . mi.,,, 


. . .  ^KSS'  ' 


102 


All  of  the  images  were  noisy  in  the  sense  that  only  a  finite  number 
of  significant  figures  were  known.  Specific  bounds  or  limits  were  not 
derived  which  would  enable  us  to  explcitly  state  the  exact  number  of 
significant  figures  used.  However,  the  number  of  significant  figures 
available  for  use  in  the  restoration  was  estimated  in  each  case. 

In  some  cases  the  image  accuracy  available  was  limited  only  by  the 
number  of  points  chosen  to  approximate  the  integration.  These  images 
were  accurate  to  at  least  4  or  more  sighificant  figures  and  are  referred  to 
as  the  "no-noise"  images.  To  obtain  images  with  less  accuracy  than  the 
no-noise  cases,  the  accurate  images  were  perturbed  by  an  additive 
random  variable. 

The  determination  of  an  estimate  for  the  available  significant 
figures  and  other  measures  of  the  amount  of  noise  used  in  the  various 
restorations  are  discussed  below. 
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a 

b(£)  =  f  sinc2(4-o)do 
-a 


which  is 


(1*9) 


b{£)  =  11 rSl  cos  2ff<fe+a)  "  (£+a)  cob  2ir(4-a) 

2ir2(£+b)(|-a) 

+  “  j  Si[2ir(4+a)]  -  Si  [2ir(£-a)]  J.  (1 70) 


where 

Si(x)  =  f‘-±Z  du. 
o 

Using  (170),  the  numerical  image  accuracy  from  the  computer  was 
compared  with  the  true  value  (170)  for  at  least  one  image  argument.  The 
minimum  number  of  significant  figures  tc  which  these  two  values  agreed 
is  defined  as  the  MSF. 

To  obtain  the  images  resulting  from  the  smooth  sources,  one  must 
perform  the  integration 

a 

b(|)  =  J' sine  (£-o)  cos  kir  (0-9)  do. 

-a 


(171) 
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This  integration  is  not  easily  performed  analytically,  nor  is  an 
estimate  of  the  error  for  Simpson's  summation  easily  obtained. 

For  these  images  a  rough  one  significant  figure  approximation  war 
determined  analytically  which  insured  that  the  image  values  were  of  the 
right  order  of  magnitude.  As  a  further  check,  the  number  of  increments 
(points)  in  the  integration  interval  was  increased  while  comparing  the 
image  values  obtained.  The  MSF  was  recorded  as  the  minimum  number 
of  significant  figures  to  which  the  two  numerical  values  agreed. 

Image  accuracy- -additive  noise  cases 

In  the  simulated  restorations,  the  noisy  image  was  obtained  using 
the  additive  model 


b  =  b  +  n  (172) 

where  n  is  a  Gaussian  random  vector  of  mean  zero  and  covariance  matrix 

2  • 
or  L  Thus,  b  is  a  Gaussian  random  vector  with  mean  b  and  covariance 
n 

2 

matrix  ar  I.  In  each  case  where  noise  was  added,  the  true  "no-noise" 
n 

image  value  (b.)  was  known  accurately  and  then  perturbed  by  the  random 

noise  sample.  Thus,  when  noise  is  not  added,  it  is  assumed  that  the  MSF 

described  previously  are  available  for  the  restoration  process. 

2 

It  is  desirable  to  determine  a  criterion  which  relates  ar  to  the 

n 

number  of  significant  figures  available  for  various  image  strengths  b.. 
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Several  criteria  could  be  established  to  find  such  a  relationship.  The 
criterion  chosen  is  explained  below  and  is  one  which  agrees  closely  with 
the  actual  simulated  results. 

Let  the  true  value  fe^  =  1.0,  and  consider  choosing  <r  ^  such  that 

P[-i<n.<l]=i  (173) 

and 

P[n.  <-iJ  +P[ui>£j  =i 

Thus 

p[i  <b.  <1.5]  =  £  {174) 

and 

P[b.<i]  +P[bl>1.5]  =  £. 

Now  consider  the  following  rounding  procedure.  (The  image  values 
were  not  actually  rounded  but  the  procedure  is  considered  here  purely 
for  establishing  a  reasonable  <r ■  ^  vs.  significant  figure  relationship. ) 
Suppose  we  set  bj  =  1.  0  if  .  5  <  b^  <  1.  5.  If  bj  <  .  5  we  set  K  =  0.  0,  and 
if  b^  >  1. 5  we  set  b^  =  2,  0,  Thus  there  is  a  probability  of  .  5  that  b.  vri*’ 


lit. 
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have  two  significant  figures  and  a  probability  of  .  5  that  will  have  no 
significant  figures.  In  this  case  the  mean  number  of  significant  figures 
is  1. 

2 

To  find  <r  we  have 
n 


•Sh*  -u2/2 
=  f  du=i 

5/<r 


(175) 


2  2 
and  <r^  =  .  55.  For  2  significant  figures  and  b.  b  1,0,  we  find  <r ^  such 


P[.  05  >  n.  >.05] 


P[n.  <  05]  +  P^  >  .  05]  =  £. 


(176) 


The  noise  variance  is  <r  =  .  055. 

n 

Figure  9  presents  log  [1  /<r  ^]  vs,  the  mean  number  of  significant 

n 

figures  for  the  various  image  strengths  used.  The  cross  hatched  region 
is  the  region  used  most  frequently  in  this  paper. 

Since  each  image  varies  in  strength  as  a  function  of  its  argument, 


then  for  <r  fixed  the  number  of  significant  figures  available  varies  over 
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2 

a  ra:.ge  of  value  This  range  corresponding  to  a  fixed  value  of  is 

defined  as  SFR,  and  when  <r  **  >  0  it  is  recorded  on  each  restoration. 

n 

As  a  further  measure  of  image  raise  level  we  use  the  following 
sample  mean  absolute  error: 


« 

8 


(177) 


To  indicate  the  relative  magnitude  of  the  noise  the  following  percentages 
are  defined, 


P 

max 


«  x  100 

s _ 

(max)  * 


ana 


a  78) 


<  x  100 

P  •  =  T -  . 

nun  6^  (min) 

In  (1 78)  tiie  sample  mean  error  is  computed  as  a  percentage  of  the  true 
"no-noise"  maximum  and  minimum  image  measurements  used  in  the 
restoration  process.  The  reciprocals  of  these  quantities  are  regarded 
as  indications  of  the  maximum  and  minimum  SNR's  used. 
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Obtaining  the  A  matrix 

In  practice  the  point  spread  function  could  be  either  measured  or 
analytically  approximated.  In  either  case  some  error  is  involved. 
However,  it  is  well  to  keep  in  mind  that  the  A  matrix  would  be  less  noisy, 
for  many  applications,  than  the  image  because  the  point  spread  function 
could  be  accurately  determined  from  many  samples  of  the  transfer 
function  in  a  controlled  low  noise  laboratory  situation. 

We  have  considered  two  general  cases  of  A  matrix  accuracy,  one 
which  uses  the  full  analytical  accuracy  of  the  computer  to  determine  the 
A  matrix,  and  the  more  impor^unt  practical  case  when  the  A  matrix  is 
noisy.  In  the  latter  case  several  levels  of  noise  were  considered  which 
exceeded  the  image  noise. 

A  matrix  accuracy 

The  additive  model  used  in  obtaining  the  noisy  A  matrix  may  be 
considered  from  two  equivalent  viewpoints.  First  suppose  we  have  a 
noisy  A  matrix  consisting  of  the  true  A  matrix  plus  a  B  matrix  composed 
of  Gaussian  random  variables.  Thus  the  image  is 

b:(A  +  E)x, 

or 

b  =  Ax  +  Bx . 


(179) 


(ISO) 


no 


Since  Bx  is  a  vector  whose  elements  are  linear  combinations  of  Gaussian 
random  variables,  the  elements  of  Bx  are  also  Gauasianly  distributed. 
Thus,  Equation  (lft0)  may  be  alternati\  -  written  as 


b  =  Ax  +  n 


(181) 


where  n  =  Bx.  This  model  is  equivalent  to  the  original  additive  model  in 
Equation  (114).  Alternatively,  if  the  A  matrix  or  the  point  spread  function 
is  determined  by  measuring  the  image  response  to  a  point  source,  then 
we  may  write  the  imaging  equatior  as 


a. .  s  a, .  +  n. 
ij  ij  x 


(18?) 


where  the  a.,  are  the  true  values  and  n.  is  the  additive  background  noise 
ij  i 

referred  to  previously. 


To  obtain  the  noisy  A  matrices  Equation  (182)  was  used  with  n.  ~ 

2  -  2  - 

N(0,  <rA  )  and  a..  -  N(a..,  c  .  ).  The  a.,  values  were  the  analytic  values 
A  ij  ij  A  ij 

2 

obtained  by  the  computer.  When  <r  "  0  the  full  analytic  accuracy  of  the 
comouter  (usually  Q  significant  figures)  was  used. 

For  the  noisy  A  matrix  cases  we  have  chosen  to  indicate  the  le/el 


of  noise  by  defining 


Ill 


N  M 
=  S  2 


i=l  j=l 


a.. -a.. 

U _ 1L 

NM 


(183) 


as  the  sample  mean  absolute  error  between  the  random  variable  a. .  and 

ij 

the  true  analytic  value  a_.  Then  we  compute 


<  A  x  100 

p  =  Jt _ . 

max.  a.,  (max)  ' 
A  ij 


(184) 


The  minimum  percentage  is  not  computed  because  the  a_(min)  values 
were  often  aero  or  very  close  to  aero.  Thus  the  effective  SNR  for  the 
noisy  A  matrix  elements  ranges  from  a  minimum  of  zero  to  a  maximum 


of  1/P 


max. 


Improvement  Criterion 


The  choice  of  a  measure  which  indicates  restoration  improvement 
over  the  initial  object  estimate  in  the  absence  of  data  processing  (which 
we  have  assumed  throughout  is  the  image)  will  likely  depend  upon  the 
application.  In  the  absence  of  a  specific  application  we  consider  two 
improvement  measures,  the  ad-hoc  or  intuitative  visual  measure  and  a 
mathematical  measure.  Generally  speaking,  these  two  measures  should 
agree,  at  least  on  a  qualitative  basis. 

In  many  problems  one  can  assess  mathematical  improvement  by 
using  the  "standard"  squared  error  (SE)  criterion  which  is 
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2 

SE  S  (x.  -  x.)  (185) 

i=l  1  1 

where  x^  is  the  estimate  resulting  from  the  restoration  process  and  x^ 

is  the  corresponding  value  of  the  true  object.  We  have  defined  the 

general  term  squared  error  (SE)  so  that  the  term  for  mean  square  error 

(MSE)  can  be  used  exclusively  as  defined  previously.  The  major  difference 

2 

between  the  MSE  and  the  SF  as  shown  above  is  that  when  <r  is  finite  we 

x 

are  inferring  the  possession  of  a  priori  information  which  we  do  not 

actually  have.  If  the  prior  guess  for  the  object  is  the  image,  then  a 
2 

finite  <r  will  weight  the  image  ao  that  the  estimate  will,  roughly  speaking, 

2 

lie  somewhere  between  the  image  and  the  true  object,  but  because  <r  is 

X 

finite  we  do  not  know  the  exact  error  involved. 

To  continue  the  discussion,  it  was  evident  that  visual  improvement 
could  not  be  assessed  only  by  SE  improvement.  In  fact,  as  will  be 
illustrated,  there  were  several  cases  when  visual  improvement  was 
apparent,  but  there  was  actually  a  decrease  in  SE  improvement.  The 
difference  is  due  to  the  fact  that  the  SE  does  not  account  for  object  size. 

The  squared  distance  (162),  on  the  other  hand,  does  include  size  error, 
and  it  agreed  more  closely  with  the  intuitative  visual  error.  Because 
of  this  agreement  we  have  chosen  to  indicate  mathematical  improvement 
by  computing  the  following  ratio 
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where  the  x.'s  anda.'s  are,  respectively,  the  true  object  values  and 

2 

arguments  and  M  a  N.  Mathematical  improvement  is  shown  if  I(d  ) 
exceeds  unity. 

Of  course  there  is  some  arbitrariness  as  to  the  choice  of  parameters 
in  (186).  For  b^,  and  M  the  values  were  taken  as  those  actually  used. 
The  x^  and  a.  were  chosen  such  that  the  were  equally  spaced  in  the  true 
nonzero  source  interval.  The  a  and  corresponding  x^  were  chosen  such 
that  the  estimated  object  was  reasonably  represented  in  the  estimated 
nonzero  source  interval.  To  clarify  this  last  statement  it  was  assumed, 
on  the  basis  of  a  posteriori  information,  that  the  regions  where  the 
estimated  object  was  negative  or  small  positively  were  regions  where  the 

estimated  object  intensity  was  zero. 

2 

The  presentation  of  I(d  )  on  each  restoration  will  allow  a  more  valid 
assessment  of  improvement  over  the  image.  Unless  otherwise  noted, 
the  I(d^)  shown  on  the  restorations  is  the  improvement  of  the  last  iteration 
over  the  image. 

It  is  well  to  remark  that  there  are  two  general  categories  of  simulates 
restorations:  those  which  assume  essentially  no  prior  size  information 
and  those  which  do  assume  some  prior  rise  information.  In  the  former 
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case  I(d  )  will  present  a  valid  improvement  criterion.  However,  for  the 

2 

latter  case  much  of  the  l(d  )  improvement  may  be  due  to  the  fact  that 
2  2 

the  term  (ar^-o^)  near  zero  while  i®  still  large.  That  is, 

2 

for  the  latter  case  it  may  well  be  that  much  of  the  X(d  )  improvement  is 
due  to  the  assumed  size  information.  Throughout  the  results  which 
follow  we  have  been  careful  to  emphasize  those  cases  in  which  prior 
size  information  is  assumed  so  that  the  reader  can  more  critically  assess 
the  restoration  improvement. 

Results  of  the  Various  Restoration  Schemes 

The  solutions  investigated  by  actual  simulations  are  identified  as 
follows: 

Solution  Description  Equation  Number  or  Reference 


MSE  (155) 

Phillips -Twomey  (44) 

Dynamic  Programming— Area  (78) 

Dynamic  Programming— Prior  Vector  (Bellman,  1965) 

Matrix  version  of  the  area  constraint  (66) 


To  identify  which  solution  was  used  and  the  various  parameters  involved, 
the  following  information  block  (Figure  10)  is  placed  on  each  restoration. 


Run 

Mi  N 

Quad  Rl  R  j  \  (  Area  {Prior  Iter 

%2 

P  P  .  ISFR  jMSF  tPmax. 

max  min  |  i  VA  r  A 

lid!) _ 

Figure  1 0.  Information  block. 
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When  the  symbol  NA  appears  in  place  of  any  of  the  quantities  in  Figure  1 0, 

it  means  that  the  quantity  does  not  apply  for  the  case  in  question. 

Run  is  an  identifying  number.  The  first  row  of  quantities  under  Run 

specifies  the  measuring  and  predicting  scheme  used  and  the  solution  used. 

The  next  row  indicates  the  image  and  A  matrix  accuracy  and  the  a  priori 
2  2 

variance  <r  when  or  is  specified.  The  last  row  indicates  the  mathe- 
x  n 

2 

matical  measure  of  improvement.  The  quantities  M,  N,  *,  R,  <r  , 

2  2 

P  .  P  .  ,  SFR,  MSF,  o- .  ,  P  ,  and  I(d  )  have  been  explained, 
max'  min  A  max. 

A 

Further  explanation  of  these  parameters,  as  applicable,  and  those  not 
previously  mentioned  is  presented  below. 

The  following  explanation  refers  to  the  parameters  Area,  Prior, 
and  Iter,  which  may  be  used  to  identify  the  type  of  solution  used.  To 
avoid  any  possible  confusion,  the  solution  description  is  also  recorded  in 
each  of  the  figure  captions.  When  Area  is  greater  than  zero  it  is  the  true 
area  tinder  the  object  and  indicates  that  the  area  constraint  wac  used. 

When  Prior  is  greater  than  zero  then  the  prior  vector  constraint  was  user*.. 
If  Prior  s2,0  then  the  image  was  used  as  the  initial  prior  vector.  When 
Prior  and  Area  are  both  zero. then  the  Phillips -Twomey  smoothing  matrix 
was  used.  The  quantity  Iter  is  the  maximum  nur.ber  ot  iterations  used 
in  the  iteration  procedure.  In  the  cases  where  Iter  >  0  and  the  Phillips- 
Twomey  smoothing  constraint  or  the  matrix,  version  of  the  area  constrair  v 


were  used,  this  is  an  indication  that  only  the  initial  prior  vector  was 
determined  by  these  methods.  Further  Iteration  used  the  MSE  solution 
(155). 

The  variance  c  ^  specifies  the  amount  of  pseudo-priori  information 

2  2  2  2 
assumed  in  the  smoothing  ratio  <r  fa  .  Figure  11  shows  cr  vs,  or 

n  X  !Q  X 

2  2 

for  various  < r  fa  ratios.  The  crosshatched  region  is  that  used  most 
n  x  6 

frequently  in  this  paper. 

The  ratio 


R1  = 


Measuring  interval  width 


Predicting  interval  width 


(187) 


relates  pertinent  information  about  the  measuring  and  predicting  scheme  r 
The  proper  choice  of  this  ratio,  as  will  be  demonstrated,  can  significant, 
alter  the  MSE  in  the  restoration  process.  In  order  to  differentiate 
between  the  cases  when  one  prediction  or  two  prediction  intervals  were 
used  for  the  same  level  of  diffraction  R,  the  prediction  interval  width 
in  the  latter  case  is  taken  as  the  sum  of  the  two  intervals. 

Three  different  quadrature  weights  were  investigated.  When 
Quad  s  1,0  the  most  extensively  used,  and  most  powerful.  Gauss 
quadrature  is  indicated;  Quad  s  2.0  corresponds  to  the  use  of  Simpson's 
quadrature,  and  Quad  :  3.0  refers  to  the  case  when  the  weights  were 


unity. 
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Effect  of  prior  information  in  the  inverse  operator 

A*  an  introduction  to  the  computational  work  a  sequence  of  figures 
is  shown  illustrating  tine  importance  of  using  finite  c r  ^  in  the  inverse 
operator. 

Figure  12  shows  the  restoration  of  a  uniform  pulse  by  operating  on  a 
low  noise  image  (MSF  *  5.  0)  with  the  straightforward  matrix  inverse  (32). 
The  esti  d  object  is  in  error  because  more  measurements  are  needed 
in  the  restoration  procedure  to  adequately  describe  the  object.  That  is. 
the  restoration  is  quadrature -error  limited. 

To  reduce  this  error  the  number  of  points  in  the  [-1,  1  ]  interval 
was  increased  from  5  to  7.  Figure  13  shows  the  results, and  the  improve¬ 
ment  is  noticeable.  Now  the  question  to  consider  in  further  improving  the 
restoration  is  the  upper  limit  of  the  number  of  points  one  can  use  in  the 
interval  [-1, 1].  The  answer  to  this  question  is  that  the  number  of  points 
for  reasonable  restoration  depends  upon  the  system  noise.  In  fact  the 

MSE  (159)  depends  upon  the  number  of  points  one  uses  and  increases 

2  . 

rapidly  as  this  number  increases.  If  the  system  noise  <r ^  is  not  small 
enough  to  overcome  this  increase  in  dimensioruJlity,  then  poor  restoration 
results.  To  illustrate  tills,  the  noise  level  of  Figure  13  was  increased. 
The  noise  variance  was  raised  to  2. 5  x  10*  which  resulted  in  a  SFR  of 
1.4-1.  7,  while  the  previous  image  accuracy  was  at  least  5  significant 
figures. 


4.  0 


pulse  (R  =  1.0)  using  (32). 
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Operation  on  the  noisy  image,  in  the  same  manner  as  before  (32), 
produced  the  oscillatory  object  estimate  shown  in  Table  1  below.  This 
estimated  object  is  meaningless.  The  imawa  noise  is  now  the  limiting 
factor.  A  similar  result  in  terms  of  an  oscillatory  solution  was  obtained 
by  increasing  the  dimensionality  to  9  and  using  the  "no  noise"  image  of 
MSF  -  5. 0. 

The  assumption  of  finite  prior  information  can  alleviate  the  oscillation. 

2  -2 

Consider  Figure  14  where  we  have  assumed  <r^  •  2.  5  x  10  and  used 
the  Phillips -Twomey  solution.  The  reduced  SE  in  Figx  re  14  over  the 
values  shown  in  Table  1  confirms  the  usefulness  of  a  priori  information 
in  the  inverse  operator. 


.300 


1.0 
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The  sequential  estimation  of  a  unifor  m  pulse 

Here  we  present  results  on  the  estimation  of  the  size  and  shape  of 

a  uniform  pulse  which  has  been  diffracted  5  times  greater  titan  the  level 

corresponding  to  the  Rayleigh  .resolution  criterion  for  two  point  sources 

(R  n  10. 0).  The  importance  of  these  figures  is  that  they  demonstrate 

the  sequential  estimation  of  an  object  which  has  been  excessively 

diffracted,  for  the  case  when  we  have  assumed  virtually  no  a  priori 

information  about  the  true  object. 

Figures  1 5-20  present  the  series.  The  initial  size  in  Figure  15 

[4. 0,  1. 0]  war  assumed  larger  than  one  would  estimate  on  the  basis  of 

comparing  widths  of  the  point  spread  function  and  image.  However, 

notice  that  the  size  is  substantially  reduced  from  [-1.0,  1. 0]  to  [  >.46, 

•  46]  •  Using  the  updated  size  estimate  the  next  figure  shows  further 

improvement  in  size  along  with  shape  improvement.  This  trend  continue  3 

until  Figure  19,  which  no  longer  indicates  noticeable  size  improvement* 

No  attempt  was  made  to  "push"  the  results  to  a  limit  by  iterating  further. 

However,  the  impending  improvement  in  >ize  and  shape  is  apparent. 

When  the  size  is  assumed  smaller  than  [  •,  2,  .2],  Figure  20  shows 

that  the  end-points  of  x  tend  toward  larger  values  than  the  true  object. 

This  suggests  that  the  actual  size  limit  of  the  restoration  process  has 

2 

been  exceeded,  at  least  for  the  particular  value  of  <r  used. 

X 


14.41 
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The  sequential  estimation  of  a  smooth  pulse 

Result*  are  presented  here  for  the  restoration  of  smooth  cosine 
pulses*  These  restorations  may  be  compared  visually  with  those 
presented  in  Figures  12-20  to  assess  the  loss  incurred  when  the  objects 
are  not  smooth  but  have  sharp  corners*  As  in  the  last  series  (Figures 
?  5-20),  these  figures  demonstrate  restoration  when  no  a  priori  object 
information  is  used* 

Figure  21  showt?  that  the  source  cos  (iror/2)  having  undergone  a 

diffraction  such  that  R  ■  1. 0  can  be  very  closely  restored.  Since  the 

diffraction  was  small  only  one  size  estimate  was  necessary. 

‘Ihe  estimation  of  the  source  cos  (2mx)  is  shown  in  Figures  22-26. 

Both  the.  diffraction  and  image  noise  are  greater  than  in  Figure  22,  but 

the  rinal  restoration  Figure  26  is  comparable  to  that  of  Figure  21.  This 
Z 

suggests  that  <tq  in  Figure  21  could  have  been  increased  with  jut 
significantly  degrading  the  restoration. 

Figures  25  and  26  are  identical  except  for  the  prior  weighting. 

Figure  26  (tr^  =  }  shows  an  increase  in  restored  detail  in  both  the 
first  and  20th  iteration*  However,  iteration  20  is  asymmetrical.  This 
feature  is  visually  noticeable  at  the  interval  endpoints  and  occurs  due 
to  computational  effects  within  the  computer.  A  computational  alternative 
is  to  choose  a  lower  value  of  <r^f  begin  with  a  more  stable  initial  solution, 
and  iterate  longer*  Here  again,  the  solution  eventually  becomes  unstable. 
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Figure  22.  Restoration  of  the  pulse  cos  (2  ir  a)  using  '155) 
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Restoration  of  the  pulse  cos  ( Z  irar  )  using  (155)  with  x 


Figure  26.  Restoration  of  the  pulse  cos  (2ir a  )  using  (155)  with  x 
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However,  we  have  the  option  of  picking  any  solution  between 
iteration  1  and  20  of  Figure  26,  Iteration  6  was  very  close  to  the  true 
source,  and  one  could  ascertain  by  visually  observing  the  plotted  estimates 
that,  as  the  number  of  iterations  increased  beyond  6,  the  computational 
inaccuracy  began  to  increase* 

The  sequential  estimation  of  split  uniform  pulses 

In  this  series  (Figures  27-32),  we  can  assess  the  importance  of  prior 

size  information.  It  is  demonstrated  thet  the  lack  of  size  knowledge  can 

seriously  affect  the  ability  of  the  restoration  process  to  restore  object 

detail.  It  will  also  be  shown  that  we  can  alleviate  this  problem  somewhat 

by  considering  various  size  estimates  and  using  an  alternate  size 

estimation  scheme  in  which  we  initially  choose  a  very  small  size  estimate 

instead  of  a  large  one  as  used  previously. 

Figures  27-28  illustrate  an  effort  to  determine  the  outside  size  of 

2  -2 

the  object.  Notice  that  iteration  20  for  <r^  =10  is  badly  distorted,  but 

for  <rx  a  10**  tiie  20th  iteration  is  stable.  These  figures  suggest  that 

the  noise  level  is  sufficiently  high  to  limit  outside  size  determination, 

2 

at  least  for  the  values  of  <r  and  the  number  of  iterations  considered. 

■x 

Figures  29-30  show  considerable  visual  improvement  over  Figures 

27-28,  This  improvement  is  attributable  to  the  assumption  of  the  true 

outside  size  in  Figures  29-30.  A  comparison  of  Figures  29  and  30  allows 

2  2  -5 

a  visual  assessment  of  the  effect  in  changing  ,  (<r^  =  10**3  for  Figure 
29  but  is  reduced  to  10~^  in  Figure  30), 


i55)  with  x 


Run 


igure  30.  Restoration 
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Figure  31  presents  a  restoration  in  which  +he  assumed  outside  size 
interval  was  smaller  than  the  true  value.  Note  the  high  values  at  the 
endpoints  of  x  which,  as  before,  suggest  that  the  size  estimate  is  too 
small.  This  figure  also  illustrates  that  the  SE  is  not  adequate  in 
measuring  improvement  because  it  is  evident  that  there  is  considerable 
gain  in  visual  information,  but  the  SE  actually  shows  an  improvement 
decrease.  These  results  add  additional  insight  to  the  size  and  shape 
estimation  procedure.  They  indicate  that  we  should,  in  any  case, 
consider  both  large  and  small  size  estimates.  Also  indicated  is  an 
alternative  size  estimation  procedure,  which  initially  chooses  a  small 
prediction  interval  close  to  the  image  maximum. 

Figure  32  shows  a  single  run  attempt  at  restoring  4  narrow  closely 
spaced  sources.  The  restoration  does  indicate  two  regions  of  enhanced 
object  intensity  which  could  be  further  investigated  by  using  split 
prediction  intervals,  as  demonstrated  below. 

The  sequential  estimation  of  two  split  smooth  pulses 

The  last  series  showed  results  of  estimating  the  outside  size  of  a 
spl*t  object  by  assuming  various  size  estimates.  Here  we  demonstrate 
(Figures  33-37)  the  use  of  split  prediction  intervals  in  estimating  the 
inside  size  of  a  split  object.  The  last  two  fi«uz  .a  (Figures  38  and  39) 
illustrate  that  we  can  resolve  two  narrow  pulses  which  have  been 
excessively  diffracted. 
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In  the  aeries.  Figures  33*37,  on  split  prediction  intervals,  Figure  33 
establishes  the  initial  inside  sice*  Figures  34-36  indicate  further 
improvement.  Notice  again  how  size  improvement  generally  implies 
shape  improvement.  Figure  37  shows  the  improvement  when  the  true 
size  is  known.  These  figures  illustrate  that,  even  though  the  objects  are 
diffracted  to  a  level  twice  as  great  as  the  Rayleigh  level,  there  may  be 
a  possibility  of  not  only  resolving  the  objects  but  of  restoring  detail  as 
well. 

a  a  follow  up  of  the  sequence  just  considered,  one  logically  seeks 
the  minimum  separation  between  two  objects  which  will  still  allow  the 
two  objects  to  be  resolved.  This  minimum  separation  was  not  found 
analytically  as  a  function  of  noise  parameters,  but  Figures  38.39  show 
results  obtained  from  the  restoration  of  two  wplit  sources  which  are 
diffracted  more  than  5  times  the  usual  Rayleigh  level. 

Figure  38  presents  the  results  for  a  relatively  high  SNR  case,  and 
Figure  39  shows  the  impending  cost  for  an  increased  amount  of  image 
noise.  Both  figures  clearly  r-w  that  the  sources  are  resolved.  Based 
on  these  results  split  prediction  intervals  can  be  used  to  further  enhance 
and  restore  the  "bright  spots. " 

Before  continuing  it  is  well  to  reiterate  that  we  have  considered  here 
a  general  restoration  procedure  which  proposes  to  restore  general  objects 
in  every  detail.  We  have  assumed  that  very  little  prior  information  is 
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available  and  have  been  careful  to  point  out  those  case*  when  prior  size 
information  was  used.  If  the  problem  were  to  find  the  minimum  aeparation 
Riven  that  two  point  sources  made  up  the  object  we  could  undoubtedly  do 
better,  in  terma  of  both  noiae  and  diffraction,  than  the  reaulta  ehown  in 
Figures  38-39. 

Restoration  using  a  noisy  A  matrix 

Results  showing  the  effect  of  various  amounts  of  A  matrix  uncertainty 
coupled  with  noisy  image  measurements  are  shown  in  Figures  40-45. 

i 

These  results  are  likely  the  most  practical  of  any  shown  thus  far  and 

i 

| 

clearly  indicate  the  possibility  of  using  the  MSE  restoration  procedure  in 
practice. 

1 

2  -s 

All  of  the  figures  have  the  same  image  noise  (<r&  =10  ,  SFR  = 

2  -5 

1.3-2. 1).  Three  A  matrix  noise  levels  are  represented:  =10  , 

2-4  2-3 

<r  =10  ,  and  <r .  =10  ,  The  first  3  figures  use  an  initial  size 

A  A 

estimate  of  [-.5,  .5]  while  the  last  3  use  an  updated  size  estimate  of 
[-.  27,  .  27]  •  Again  we  note  the  marked  general  improvement  in  shape 
when  the  size  estimate  is  closer  to  the  true  size. 

The  first  two  figures,  Figures  40  and  41,  have  the  same  A  matrix 
and  image  noise  but  illustrate  the  computational  advantage  of  using  a 

5 

smaller  <r  ^  (<r  ^  «  1 0~^  in  Figure  41  while  r  ^  =  lO*"5  in  Figure  40)  to 

X  X  X 

obtain  a  more  stable  initial  solution,  which  can  be  iterated  longer  and 
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thus  produce  a  better  restoration.  Figure  42  then  shows  the  cost  of 

increased  noise  on  the  improved  solution  of  Figure  41. 

Figures  43-44  show  the  same  cost  as  Figures  41-42  when  the  noise 

2  .5  2  *4 

is  increased  frome^  «  10  to  <r ^  a  10  ,  except  that  the  size  is  now 

assumed  to  be  closer  to  the  true  size.  The  improvement  is  noticeable 

in  Figures  43-44  over  Figures  41-42. 

The  last  figure  in  this  series,  Figure  45,  illustrates  the  further 

cost  in  restoration  detail  when  the  A  matrix  noise  is  further  increased 

to  V  a  10  . 

A 

Image  smoothing  for  the  low  SNR  Case 

Now  we  consider  the  question  of  what  happens  when  we  encounter 
excessive  image  noise.  Suppose  the  image  we  ha*  j  at  our  disposal 
appears  to  be  too  noisy  to  improve.  We  consider  here  a  possible  approach 
to  this  problem  by  smoothing  the  original  noisy  image  in  accordance  with 
prior  information  inherently  available  due  to  diffraction  effects. 

Figure  46  presents  a  restoration  resulting  from  the  operation  on  a 
very  noisy  image  (SFR  =  .  6-1. 6).  Iteration  1  is  not  stable  for  <r  ^  - 
2.  5  x  10  In  fact  we  must  weight  the  prior  noisy  image  so  heavily  in 
order  to  stabilize  the  solution  that  there  is  essentially  no  improvement. 
Apparently  the  noise  level  f  i  so  high  that  restoration  efforts  are  futile. 

Of  course,  if  additional  samples  of  the  image  were  available,  some 
improvement  may  be  possible. 
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On  examining  “  •  noisy  image  it  is  obvious  that  the  spatial  extent  of 
the  noise  fine  structure  is  small  compared  with  the  width  of  the  point 
spread  function.  This  suggests  that  the  noise  fine  structure  in  this  case 
can  be  smoothed,  since  it  could  not  logically  be  preserved  by  the  imaging 
system.  A  resulting  smoothed  image  is  shown  in  Figure  47.  (No 
sophistication  was  used  in  the  smoothing  process. ) 

Figure  48  shows  the  results  of  using  a  9x9  system  and  operating  on 
tiie  smoothed  image.  The  improvement  is  noticeable.  These  results 
illustrate  that  even  in  a  low  SNR  situation  in  which  the  detector  noise  is 
appreciable  some  measures  can  be  taken  which  may  result  in  restoration 
improvement. 

Quadrature  effects 

As  previously  discussed  in  the  Imaging  Equation  section,  the  transfer 

from  the  continuous  functions  to  the  discrete  version  involves  some  error. 

The  visual  effect  of  this  error  is  illustrated  in  Figures  49-51.  The 

restoration  shown  in  Figure  49  used  unity  weights,  Simpson's  quadrature 

was  used  in  Figure  50,  and  Gauss  quadrature  was  used  in  Figure  51* 

Under  the  stress  of  both  image  and  A  matrix  noise  it  is  evident  that 

the  more  powerful  Gauss  quadrature  is  superior.  There  may  be  some 

question  about  the  small  values  of  x  in  Figure  49;  however,  it  was  found 
2 

that  decreasing  <r^  ,  in  general,  has  a  greater  scale  auect  for  =  1.0 
than  when  one  of  the  other  quadratures  was  used.  (This  scale  change  is 
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further  evident  if  we  consider  the  problem  of  numerically  integrating  to 
find  the  area  under  a  uniform  pulse  of  unity  height  while  using  unity 
weights*) 

The  choice  of  quadrature  and  its  relationship  to  the  MSB  is  discussed 
later  in  the  section  entitled  *'MSE  Variation— Choosing  the  A  Matrix 
Parameters* " 

Comparison  of  the  Phillips  -Twomey  and  MSE  Solutions 

If  we  accept  the  use  of  successive  approximations,  then  the  essential 
difference  between  the  Phillips -Twomey  solution  and  the  MSB  solution  is 

■] 

just  the  initial  prior  vector  used.  Figure  52  shows  the  asymptotic 
convergence  of  the  SE  vs.  iterations  for  the  two  solutions.  (The  SE  is 
increasing  with  increasing  iterations  since  thr  5x5  system  used  is 
approaching  the  quadrature  error  limited  solution  for  ^  =  0  shown  in 
Figure  12.) 

I 

The  figure  indicates  that  the  Phillips -Twomey  solution  has  a  lower 

i 

SE  than  the  MSE  solution  for  iteration  1 ;  however,  after  2  iterations  the 
two  solutions  differ  only  slightly  in  the  rate  of  convergence.  In  fact,  it 
is  merely  a  matter  of  definition  as  to  which  solution  reaches  the  asymp¬ 
totic  SE  first.  That  is,  iteration  1  for  the  MSE  solution  could  have  been 
defined  as  the  zeroith  solution,  thus  shifting  the  entire  MSE  solution  curve 
1  iteration  to  the  left. 
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Dynamic  programming,  restorations 

Figures  53-56  illustrate  the  use  of  dynamic  programming  in  restoring 
a  uniform  pulse  source  when  R  =  1 . 0.  The  prior  vector  solution  (Bellman, 
1965)  was  used  for  Figures  53-55.  Notice  that  the  solution  for  small 
prior  weighting  (Figure  53)  is  not  symmetrical  and  is  poorly  behaved  as 
one  progresses  toward  the  initial  estimates  on  the  left.  This  results 

because  the  entire  A  matrix  is  not  used  until  the  last  scalar  estimate  x.  . 

hi 

is  found. 

The  solutions  for  the  problem  Bellman  (1965)  considered  appear  to 
be  more  regular  for  the  initial  x.'s  than  the  above  results  shew.  However, 
the  prior  vector  he  assumed  was  closer  to  the  true  solution  for  the  initial 
estimates.  This  would  tend  to  reduce  the  initial  irregular  behavior. 

Larger  values  of  prior  weighting  improve  the  solution,  as  shown  in 
Figures  54  and  55,  and  it  appears  that  a  value  of  x  between  1  and  1 0  would 
result  in  even  greater  improvement. 

Results  of  the  area  constraint  dynamic  programming  solution  are 
depicted  in  Figure  56.  These  results,  along  with  results  from  the  matrix 
inverse  solution  using  the  area  constraint,  indicate  that  the  area  constraint 
is  less  stringent  in  controlling  the  characteristic  oscillatory  solutions 
even  when  prior  weighting  is  used. 

Based  on  the  results  obtained,  it  appears  that  the  dynamic  program, 
ming  solutions  investigated  are  less  effective  for  optical  restoration  in 
one  dimension  than  the  MSE  and  Phillips -Twomey  solutions. 
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Figure  53.  Restoration  using  the  dynamic  programming  version  of  (47),  (Bellman,  1965) 
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MSE  VARIATION— CHOOSING  THE  A  MATRIX  PARAMETERS 


This  section  considers  the  study  of  the  MSE  and  the  insight  the  MSE 
variation  provides  in  enabling  sound  computational  procedures  to  be 
determined.  We  first  present  further  discussion  which  defines  the  A 
matrix  parameters  for  the  computational  scheme  we  wish  to  consider. 

Next  we  present  and  discuss  the  MSE  variation  with  and  without  a  priori 
information.  The  final  topic  presents  the  procedures  which  have  been 
developed  and  used  in  the  preceding  restorations  for  choosing  the  A  matrix 
parameters. 

The  form  of  the  MSE  we  consider  here  is  that  found  by  letting 
Z  1 

K  =cr  I  and  K  =  tr  I  in  Equation  (152).  The  resulting  MSE  is 

A  A  U  IX 


MSE  «  <r 
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2  ?  *1 ' 
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(180) 


When  the  noise  variance  «r^  is  fixed,  which  is  usually  the  case  since  it 

is  determined  by  the  experiment,  there  are  essentially  two  ways  to  further 

2 

minimize  the  MSE.  We  can  either  reduce  <r  (assume  more  pseudo- 

x 

prior  information)  or  adjust  the  parameters  that  govern  the  A  matrix. 

Obviously  one  seeks  the  best  measuring  and  predicting  scheme  which 
would  entail  optimal  choices  of  M  and  N,  the  £.*s,  the  o^*s  and  the 
quadrature.  This  is  indeed  a  difficult  problem.  In  fact  such  optimal 
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choices  would  appear,  intuitively,  to  depend  upon  the  assumption  of  more 
prior  object  information  (for  example,  concentrating  on  a  certain  class 
of  well  defined  objects)*  In  this  paper  we  do  not  consider  that  such 
a  priori  information  is  available.  Here  we  have  sought  to  establish 
procedures  for  the  general  case.  To  accomplish  this  we  have  numerically 
studied  the  MSE  as  the  various  parameters  were  varied. 

A  Matrix  Parameters 

There  is  an  uncountably  infinite  number  of  one -dimensional  schemes- 
we  could  consider.  \7e  have  most  extensively  studied  the  scheme  which 
uses  M  s  N  equally  spaced  measurement  locations  (£.*8)  and  M  =  N  equally 
spaced  prediction  points  (a  7  s).  The  case  when  unity  quadrature  weights 
were  v  jed  was  most  extensively  studied,  but  results  have  also  been 
obtained  for  Simpson's  and  Gauss  quadratures.  (The  ar.'s  were  spaced 
in  accordance  with  the  Gauss  method  when  Gauss  quadrature  was  used. ) 

The  ratios  R  and  R1  were  used  to  relate  the  diffraction  level  to  the 
measuring  and  predicting  intervals  being  considered.  In  this  case  R  is 
defined  as  the  ratio  of  the  point  spread  width  to  the  prediction  interval 
being  considered. 

We  note  then  that  the  A  matrix  parameters  are  R,  Rl,  Ivl  and  the 
quadrature  weights.  Specification  of  these  parameters  determines  the 


A  matrix. 


177 


MSE  vs»  A  Matrix  Parameters  for  Largo  A  Priori  Uncertainty 

Figures  57-59  show  the  MSE  vs.  R  relationship  for  various  values 

of  R1  when  o'  ^aoo.  and  w.  *  1.0,  i  «  1.....M.  The  MSE  was  normalized 

with  respect  to  noise  in  Figures  57-58  by  considering  MSE/<r^.  The 

MSE  in  Figure  59  was  normalized  with  respect  to  dimensionality  as  well 

2 

as  noise  by  considering  MSE/<r  M. 

n 

Notice  that  the  MSE  is  lowest  when  we  measure  "over"  the  prediction 
points  (R1  a  1.0)  for  small  values  of  diffraction  (the  cases  Austin  (1966) 
treated  when  A  a  I),  but  when  R  >  3. 8  then  R1  a  2.  0  is  "best. "  Barnes 
(1966b)  was  able  to  prove  that  the  A  matrix  is  positive  definite  when  R  is 
finite  and  R1  =  1,0,  which  demonstrates  that  the  £.'s  and  cr.'s  for  any 
N  a  M  can  be  chosen  such  that  the  MSE  is  finite.  However,  as  the  above 
results  indicate,  this  does  not  preclude  the  choice  of  R1  other  than  unity 
to  further  reduce  the  MSE. 

Another  important  and  obvious  feature  can  be  deduced  by  noting  the 
marked  MSE  increase  for  the  curves  of  Figure  58  when  compared  with  the 
curves  of  Figure  57.  Figure  59  shows  the  absolute  MSE  increase  for 
R1  a  1,0  when  M  increases  from  2  to  3.  This  increase  is  most  apparent 
for  R  >  1.  0  and  is  attributable  only  tJ  the  addition  of  one  more  row  an'1 
column  in  the  A  matrix.  Y/hen  we  extrapolate  to  matrices  of  larger 
dimension  (say  24  x  24  as  used  several  times  in  the  simulations)  the  MSE 
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for  R  >  1  and  the  values  of  R1  shown  becomes  astrcnomical.  This  result 

again  empht.'ires  the  amplifying  effect  that  excessive  diffraction  3  on 

the  initial  arbitrary  noise  1  vel  and  essentially  explains  the  reason  for 

the  oscillatory  solutions  when  the  A  matrix  dimensionality  is  large  and 
2 

c r  is  infinite.  The  next  figures  we  present,  which  are  explained  below, 

X 

2 

illustrate  how  finite  <r  alleviates  this  effect. 

x 

MSE  vs.  A  Matrix  Parameter e  for  Finite  /  priori  Information 

Figure  60  illustrates  the  MSE  behavior  for  various  parameters  as 

the  prior  information  is  varied.  The  general  behavior  agrees  with 

2 

intuition  since  as  prior  information  increases  decreases)  we  infer 
the  possession  of  more  knowledge  about  the  object  and  the  MSE  decreases 
accordingly. 

Figure  61  shows  the  measuring  and  predicting  intervals  for  the  top 
four  curves  shown  in  Figure  60.  After  viewing  Figure  61  refer  again  to 
Figure  60.  Now  consider  the  asymptotic  MSE  for  the  top  four  curves. 
These  top  four  curves  show  the  advantage  of  choosing  a  measuring  interval 
which  is  wider  than  the  prediction  interval  (R1  is  large)  since  the  MSE 

7 

decreases  by  at  least  a  factor  of  10  while  the  diffraction  level,  weights, 
and  M  remain  fixed.  The  bottom  two  curves  show  a  similar  reduction  for 
a  smaller  diffraction  level. 
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sine  (u) 


Rl=.  5 


Rl*l.  0 


R1  =2.  0 


Rl  =  5.  0 


Figure  6l.  Measuring  and  predicting  intervals  superimposed  with 
the  point  spread  function  for  the  top  four  curves  in  Figure  60.  The 
symbols  MI  and  PI  represent,  respectively,  the  measuring  and 
predicting  intervals. 
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Figure  60  exhibits  results  for  the  hi  *  5  case.  Figure  62  presents 
curves,  in  the  same  manner  as  Figure  60,  for  the  simulated  results  shown 
in  the  previous  section.  From  Figure  62  it  is  evident  that  the  prior 
weighting  necessary  to  produce  the  simulations  was  great  enough  so  that 
we  are  operating  in  the  region  where  the  curves  blend  together,  as  shown 
in  Figure  60.  In  this  region  choosing  III  large  still  has  an  effect  on  the 
1-/1SE,  but  it  is  not  nearly  as  pronounced  as  for  the  asymptotic  region  of 
the  curves.  However,  the  computational  error  involved  in  the  matrix 
inversion  process  was  noticeably  reduced  by  choosing  R1  large. 

Choosing  the  A  Matrix  Parameters 

Now  that  we  have  presented  curves  showing  the  general  behavior  of 

the  MSE  with  the  A  matrix  parameters  we  discuss  the  guidelines  for 

choosing  these  parameters.  These  choices  depend  not  only  upon  the  iviSE 

but  upon  other  errors,  the  most  prominent  being  the  quadrature  error, 

2 

the  computational  error,  and  the  error  caused  by  assuming  <r^  too  small. 
In  choosing  these  parameters  it  is  evident  that  we  must  malce  these  choices 
regardless  of  the  diffraction  and  noise.  Thus,  although  we  seek  to  improve 
the  distorted  object  in  the  presence  of  these  effects,  the  diffraction  and 
noise  must  be  regarded  as  fixed  quantities  when  making  these  choices. 
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Figure  62.  Curves  of  MSE/cr  ^  vs,  the  ratio  cr  ^ /cr  ,  showing  the 

n  n  x 

region  used  in  the  simulated  restorations. 
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Choosing  the  quadrature 

The  determination  of  which  quadrature  is  "best"  or  better  than 
another  is  a  difficult  problem.  One  cannot  accurately  estimate  the 
quadrature  error  without  knowledge  of  the  integrand,  which  for  the 
restoration  problem  we  have  considered  is  unknown. 

Before  discussing  the  quadrature  choice  we  briefly  review  the 
quadratures  we  have  considered.  For  a  more  comprehensive  treatment 
the  reader  is  referred  to  Kopal  (1961),  Krylov  (1962)  or  the  more  tutorial 
McCracken  and  Dorn  (1964). 

The  Gauss  quadrature  method  is  a  very  powerful  tool  to  use  in 

numerical  integration.  Here  we  refer  to  power  as  describing  the  number 

of  points  necessary  in  the  approximating  summation  to  provide  a  specified 

accuracy.  To  illustrate  just  how  powerful  the  Gauss  method  is,  consider 

that  a  polynomial  of  degree  m  approximates  the  integrand  to  a  specified 

accuracy.  Now  suppose  we  use  the  Gauss  system  with  a  summation  limit 

of  M  points.  The  Gauss  system  is  sufficiently  powerful  that  when 

m  *t  2M-1  the  integration  is  performed  exactly  with  no  error.  For  example, 

suppose  we  decide  that  the  integrand  is  sufficiently  approximated  by  a 

polynomial  of  degree  49.  This  means  that  only  25  points  are  required  to 

perform  the  integration  exactly.  Simpsons  summation,  on  the  other  hand, 

rd 

is  only  capable  of  integrating  3  degree  integrands  exactly.  Herein  lies 
the  major  advantage  of  Gauss  quadrature  over  Simpson1  s  quadrature  (and, 
as  the  above  references  point  out,  over  essentially  all  other  quadrature  ;), 
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The  only  disadvantage  of  Gauss  quadrature  in  comparison  with  other 
quadratures  is  one  of  complexity*  The  complexity  is  due  to  the  use  of 
unequally  spaced  integrand  arguments.  This  disadvantage  is  not  very 
significant  in  the  optics  problem  because  we  can  predict  the  object  for 
unequally  spaced  points  just  as  well  as  equally  spaced  points.  The  only 
possible  drawback  would  be  the  specification  of  the  point  spread  function 
at  these  points.  If  the  point  spread  function  is  analytically  approximated, 
this  would  not  be  a  problem.  For  the  case  when  the  point  spread  function 
is  measured  some  extra  care  may  be  necessary  to  ensure  that  the  point 
spread  function  is  accurately  specified  for  all  ranges  of  its  argument. 

Now  that  we  have  briefly  reviewed  the  quadratures  considered  we 
continue  on  with  the  discussion  of  quadrature  choice. 

Recall  that  the  total  error  is  comprised  of  both  MSE  and  quadrature - 
error.  These  two  errors  are  both  related  to  the  system  dimensionality. 
We  have  seen  that  increased  dimensionality  implies  increased  MSE,  but 
it  is  well  known  that  increased  dimensionality  reduces  quadrature-error. 
Thus,  it  is  logical  to  asstime  that  a  choice  of  quadrature  should  involve 
the  trade-off  between  MSE  and  quadrature -error.  The  following  results 
illustrate  this  trade-off. 

2 

Refer  again  to  Figures  49,  50,  and  51.  The  MSE/tr^  for  the  Gauss 

2  2 
case  is  1.43  x  10  ,  for  the  Simpson’s  case  1.  26  x  10  ,  and  for  unity 

2 

weights  1.23  x  10  .  The  fact  that  the  Gauss  quadrature  resulted  in  a 
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greater  MSE/<r  2  is  easily  explained.  In  all  three  cases  the  actual  number 
n 

of  points  used  was  equal,  but  the  effective  dimensionality  of  the  Gaues 
system  is  much  greate?  than  either  of  the  other  quadratures.  Similarly, 
the  effective  dimensionality  in  the  Simpson's  case  is  greater  than  the  unity 
weights  case.  Visually  (which  is  a  good  measure  of  total  error)  the  Gauss 
system  produces  a  better  restoration.  Since  the  MSE  actually  increased, 
but  only  :  tly,  the  discrepancy  between  visual  error  and  MSE  can  . 

logically  be  attributed  to  a  reduction  in  quadrature  error.  In  fact,  it 
appears  that  the  reduction  in  quadrature  error  more  than  compensates 
for  the  slight  increase  in  MSE.  Thus,  in  this  case  one  would  choose 
Gauss  quadrature.  Furthermore,  we  can  infer  from  Figures  60  and  62  and 

I 

other  numerical  evidence  that  as  long  as  one  uses  sufficient  a  priori 
information  to  be  in  the  region  where  the  curves  blend  together,  use  of 
the  Gauss  quadrature  will  consistently  more  than  compensate  for  the  slight 
increase  in  MSE.  And  in  this  region  it  is  recommended  that  Gauss 
quadrature  be  used. 

On  die  other  hand,  as  v  *  approaches  infinity  die  increased  effective 

X 

dimensionality  for  the  more  powerful  quadratures  may  increase  the  MSE 
above  the  compensating  reduction  in  quadrature  error.  For  this  region 
of  small  a  priori  information  we  could  either  reduce  the  Gauss  system 
dimensionality  or  perhaps  successfully  use  a  less  powerful  quadrature. 


189 


Choosing  toe  ratio  of  the  measuring  and  predicting  intervals  (Rl| 

The  infraction  range  considered  in  this  paper  is  R  >  1«  0.  For 
R  <  1, 0,  the  image  becomes  a  much  better  prior  guess  for  die  object, 
and  we  can  heavily  weight  the  image  in  the  restoration  procedure. 

Consider  the  images  for  the  Range  R  >  1.  0.  For  the  noise  levels 
considered  in  this  paper  the  image  width  (noise  to  noise)  is,  grossly 
speaking,  the  width  of  the  point  spread  function.  Thus,  the  prediction 
a  id  measuring  intervals  will  be  "under"  the  main  lobe  of  the  point  spread 
function.  When  this  is  the  case,  as  has  been  demonstrated  in  Figures  57, 

5C,  and  60,  it  is  profitable  to  choose  R1  larger  than  unity  so  that  the 
measuring  interval  is  wider  than  the  prediction  interval.  The  following 
discussion  illustrates  what  is  happening  to  the  A  matrix  elements  as  R1 
changes  and  then  presents  a  guide  for  choosing  an  upper  bound  for  Rl. 

Refer  again  to  Figures  57  and  58.  The  tendency  for  A  to  approach  a 
singular  condition  for  Rl  =  1.0  and  increasing  R  is  clearly  shown.  As 
one  studies  the  A  matrix  elements  for  the  values  of  R  >  1.  0  and  Rl  *  1.0, 
it  is  seen  that  they  become  more  and  more  alike.  On  the  other  hand, 
choosing  Rl  >  1.  C  for  values  of  R  >  1.  0  tends  to  make  the  A  matrix  elements 
differ,  and  generally  speaking,  if  Rl  is  not  made  too  large,  the  MSE  will 
decrease.  The  choice  of  an  upper  bound  for  Rl  involves  a  trade-off 


between  two  factor  3 
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The  first  is  seen  in  Figure  57.  For  matrices  of  small  dimension 
it  is  evident  that  the  A  matrix  will  be  singular  for  certain  choices  of  R1 
even  when  R  >1,0.  However,  it  is  evident  intuitively  that  if  the  dimension¬ 
ality  is  large  enough  (3x3  or  greater)  and  R1  is  chosen  only  large  enough 
so  that  the  measuring  and  predicting  intervals  both  lie  within  the  point 
spread  function  width,  then  the  A  matrix  will  not  be  singular.  An  analytic 
proof  of  this  statement  is  not  given,  but  the  available  numerical  evidence 
shows  that  A  is  nonsingular  for  this  choice  of  Rl. 

The  second  limitation  on  large  Rl  is  imposed  by  system  noise  (or^). 
For  R  >  1,  0  it  is  evident  that  as  Rl  increases  above  unity  then  the  image 
is  being  measured  where  the  SNR  is  decreasing.  Thus,  if  the  system  noise 
is  great  enough  a  choice  of  Rl  extended  to  the  limit  described  above  may 
be  suboptiroum. 

Under  these  conditions  the  following  guide  for  choosing  Rl  is  stated. 

For  matrices  of  dimensionality  3x3  or  greater  and  R  >  1. 0  choose  Rl  as 

large  as  possible,  while  ensuring  that  the  measuring  interval  is  "under" 

the  point  spread  function  and  small  enough  that  the  measurements  are  not 

excessively  noisy.  This  notion  is  clearly  a  compromise  between  system 
2 

noise  (<r  )  or  SNR  and  the  location  of  ♦he  measurements  such  that  the 

n 

matrix  (A*  A  +  <r  ^ h  ^1)  is  non-singular  enough  to  be  successfully 
n  x 

inverted. 


a 


191 


Choosing  A  matrix  dimensionality 

Throughout  the  paper  we  have  assumed  that  the  A  matrix  is  of 
dimensionality  NxM  and  that  M  <  N.  Y/h  a  M  <  N  we  have  the  over¬ 
determined  case  in  which  we  measure  the  image  more  often  the 
number  of  unknowns  we  have  to  predict.  In  some  problems  it  is 
advantageous,  for  statistical  reasons,  to  choose  M  <  N  and  use  the  over- 
determined  solution.  This  case  was  considered  in  the  simulation  presented 
in  Figure  46.  Recall  that  the  image  was  excessively  noisy  (SFR  a  .6-1.6). 
In  an  effort  to  restore  the  object  several  A  matrices  with  M  <  N  were  used. 
The  dimensionality  shown  in  Figure  46  was  50x11.  However,  in  all  of  the 
cases  when  M  <  N  no  appreciable  improvement  was  apparent.  Based  on 
this  evidence,  it  appears  that  choosing  M  <  N  is  of  no  particular  advantage 
in  the  restoration  problem.  When  M  =  N,  then  the  choice  of  how  large 
to  make  N  or  M  is  to  be  considered.  This  choice  is  discussed  as  follows. 

As  previously  stated,  the  MSE  is  proportional  to  N,  and  the  quadrature 
error  is  inversely  proportional  to  N.  Thus,  N  should  be  chosen  by 
studying  the  trade-off  between  MSE  and  quadrature  error.  For  the  one¬ 
dimensional  restorations  considered,  «r  ^  was  assumed  small  enough  so 
that  reasonably  large  matrices  could  be  successfully  uoed.  The  largest 
matrix  used  was  a  65x65,  and  it  resulted  in  restoration  Improvement. 

For  tiie  one -dimensional  objects  considered,  it  appeared  that  the  48x48 
Gauss  system  was  adequate  for  reducing  quadrature  error. 
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To  obtain  the  same  quadrature  error  for  two  dimensional  objects,  the 
A  matrix  dimensionality  will  increase  roughly  as  N2  compared  to  N.  Thus 
for  two  dimensional  objects  we  will  have  to  consider  more  closely  the 
trade-off  between  MSF  and  quadrature  error. 


2 

Choosing  <r^  and  the  number  of  iterations 

2 

The  choice  of  <r^  involves  the  compromise  between  MSE  and  the 

error  involved  by  assuming  too  much  prior  information.  One  desires  to 
2 

choose  small  enough  to  enable  stable  computational  results  to  be 
obtained  and  yet  large  enough  to  reduce  die  smoothing  effect  or  a  priori 
image  weighting. 

2 

The  simulated  results  have  indicated  that  <r^  may  be  sequentially 

2 

determined  by  visually  judging  restorations  for  several  values  o*  ^  . 

Xt  has  also  been  demonstrated  that  one  can  circumvent  the  need  v 

2 

continually  invert  matrices  for  each  <r^  in  order  to  vary  prior  information. 

Prior  information  can  be  easily  varied  by  using  successive  approximations, 

Considering  that  the  other  parameters  are  fixed,  a  good  procedure  to 
2 

use  is  to  vary  <r  until  the  solution  appears  (visually)  to  be  stable  and 
then  iterate  until  computational  error  limits  the  restoration. 
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SUMMARY  AND  CONCLUSIONS 
Synopsis  of  the  Paper 

In  this  section  a  brief  synopsis  of  the  entire  paper  is  presented  and 
the  important  results  are  stated. 

The  basic  problem  considered  is  one  of  restoring  an  optical  object 
which  has  been  diffracted  and  corrupted  by  noise.  Previous  work  has 
shown  that  if  the  object  is  of  finite  extent  then  the  restoration  process 
is  limited  by  noise  and  not  diffraction.  After  reviewing  previous 
approaches  to  the  problem,  it  was  evident  that  improvement  could  be 
made  by  providing  for  noise  in  toe  restoration  procedure. 

In  order  to  present  some  insight  into  toe  diffraction  process  toe 
imaging  equation  was  derived.  The  form  of  this  equation  indicates  that 
a  linear  transformation  of  the  object  constitutes  toe  basic  image -object 
relationship  and  that  diffraction  implies  the  obscuration  of  object  detail. 
Since  numerical  techniques  were  to  be  u ,ed,  the  discrete  version  of  toe 
imaging  equation  was  presented  and  toe  quadrature  error  was  introduced. 
Next  toe  straightforward  no-noise  matrix  inverse  solution  to  the  problem 
was  presented.  Two  major  difficulties  in  using  this  solution  are  that 
excessive  image  accurar  /  may  actually  be  necessary  to  effect  the 
solution  and  that  for  the  diffraction  range  of  interest  (R  >  1. 0)  toe  A 
matrix  is  nearly  singular  and  is  difficult  to  invert. 


The  next  section  presented  previous  noiseless  work  by  Phillips  (1962), 
Twomey  (1963)  and  Bellman,  et  al,  (1964,  1965)*  Their  work  demonstrated 
that  a  priori  information  in  the  form  of  a  constraint  could  be  successfully 
used  to  alleviate  the  difficulties  just  mentioned. 

In  order  to  introduce  and  discuss  the  uncertainty  caused  by  noise, 
analytical  results  for  a  detection  example  and  for  the  error  variance 
obtained  in  estimating  tho  separation  between  two  point  sources  were 
presented.  Here  the  basic  trade-off  between  diffraction  and  noise  in 
regaining  obscured  object  information  was  presented  and  discussed  in 
detail.  The  basic  trade-off  is  that  diffraction  essentially  amplifies  the 
prevailing  noise  level,  and  in  order  to  regain  object  information  for 
excessive  diffraction  the  SNR  must  be  increased. 

Turning  again  to  the  restoration  problem,  object  estimation  in  the 
presence  of  additive  and  detector  noise  was  discussed.  Both  the  Baysean 
and  MSE  approaches  in  the  estimation  of  the  object  were  considered.  In 
using  tiie  Bayes1  approach  it  was  necessary,  because  of  mathematical 
tractability,  to  use  the  suboptimum  scheme  of  operating  on  an  image 
estimate  to  perform  the  restoration.  The  Baysean  approach  did  provide 
insight  as  to  how  a  priori  information  enters  into  the  solution,  but  the 
euboptimal  scheme  did  not  indicate  the  use  of  a  priori  information  in  the 
inverse  operator.  On  the  other  hand,  the  MSE  approach  did  indicate 
how  a  priori  information  can  be  used  in  the  inverse  operator  and 
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furthermore  indicated  that  the  smoothing  process  in  the  noiseless  case 
treated  by  Phillips  (1962)  is  actually  a  trade-off  between  the  noise  level 
and  a  priori  information*  Also  discussed  in  the  Object  Estimation  section 
was  the  assumption  of  high  SNR*  Using  this  assumption  the  Baysean 
estimates  for  Gaussian  statistics  are  equivalent  to  the  MSE  estimate. 

In  the  next  section  simulated  object  restorations  were  presented. 

Here  various  parameters  were  defined  which  indicated  how  well  the 
restoration  process  performed  for  various  objects  and  for  various 
diffraction  and  noise  levels*  These  results  demonstrated  that  object 
restoration  is  possible*  The  amount  of  restoration  improvement  is 
dependent  upon  both  diffraction  and  noise.  A  specific  statement  is  that, 
generally  speaking,  improvement  is  possible  for  the  diffraction  range 
1.0  < R  <10  when  the  sample  mean  noise  level  is  as  high  as  1  to  3 
percent  of  the  image  maximum. 


Next,  numerical  results  on  the  variation  of  the  MSE  were  presented 

2 

and  discussed.  For  fixed  e*^  and  diffraction  in  general  the  only  recourse 
to  further  reduce  the  MSE  is  to  consider  varying  the  A  matrix  parameters 
or  the  a  priori  variance  <r  2.  The  following  guidelines  were  evident  from 
the  results  presented  in  this  section.  The  quadrature  used  most  profitably 
was  the  Gaussian  system.  For  the  diffraction  range  considered  the 
measuring  interval  should  be  chos - oer  equal  to  or  greater  than  the 


predicting  interval.  The  upper  limit  of  the  measuring  interval  width 
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M 


should  be  chosen  such  that  it  does  not  exceed  the  point  spread  function 

width  or  include  excessively  noisy  image  measurements.  It  is  advisable 

to  use  M  *  N.  The  choice  of  an  upper  1-mit  cn  N  involves  the  trade-off 

between  the  error  caused  by  excessive  weighting  of  the  a  priori  image 

2 

and  the  improvement  in  quadrature  error.  The  a  priori  variance  <r^ 
should  be  chosen  small  enough  to  obtain  an  initial  stable  solution  which 
can  be  successively  iterated  until  comp  ational  accuracy  limits  the 
restoration.  These  results  when  applied  to  the  MSE  estimate  and  coupled 
with  the  sequential  size  estimation  procedure  discussed  previously 
constitute  the  basic  computatic  *  scheme  dev^oped  *  »  his  paper. 


Future  Research 


The  outstanding  applied  research  need  is  to  actually  perform  a 
restoration  experiment  which  uses  the  procedures  developed  in  the  paper. 
Such  an  experiment  would  provide  information  on  the  more  real  noise 
levels  encountered  in  practice  and  should  lead  to  the  intelligent  use  of 
repeated  samples  of  the  noisy  image  and  A  matrix. 

The  general  theoretical  or  analytical  solution  to  the  problem  remains 
to  be  obtained.  To  solve  this  problem  it  appears  that  we  need  to  specify 
eigenfunctions  and  eigenvalues  for  a  general  inverse  imaging  kernel,  as 
discussed  by  Barnes  (1966a).  This  is  indeed  a  difficult  problem,  but 
certainly  object  restoration  will  not  be  complete  until  this  information  is 
available. 
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Throughout  the  simulated  results  we  have  considered  that  the  noise  is 
white  and  Gaussian*  Band  limited  noise  should  also  be  considered*  When 
bandlimited  noise  is  considered*  perhaps  the  diffraction  effects  will  not 
be  as  severe  as  those  shown  in  this  paper*  However*  it  should  be 
mentioned  that  the  computational  error  is  not  bandlimited* 
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APPENDIX 

* 

Two-Dimensional  Ga.ua a  Quadrature 


In  two  dimensions  the  imaging  equation  is 


Vi 


j'  (a,P)dflfdP 


*2*1 


The  two-dimensional  Gauss  quadrature  tor  (189)  can  be  written 
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Pm  =  *(V12)+*<Va2h'm 

and  the  quantities  y  ,  y  ,  H  and  H  are  defined  from  tables  of 

n  m  m  n 

coefficients  for  a  given  N  (Krylov,  1962). 
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(192) 


Y'e  define 


*vi 


N  M 


tM<b,«).mtoRMW.n^&l(riWlx1<)  +t*«  *  Vj-ty  )•  (193) 

•  • 

j{  j{ 

1  X1 

According  to  the  optimality  principle  (stepwise  optimization  principle)  we 
can  write, 


M 


2  N  hi  2 

£-.(!>, c)»min[min h(S  wx-c)  +  S  (S  a  x  -b )  ] 

M  XM-1  i»l  i=l  j=l  3  3 


(194) 


and  theze  e:  1st  arguments  arg^  and  such  that  the  term  in  brackets 
can  be  written 


r  M-l  2  N  M-l  f) 

W“*l*  "*2>  *  ,»  +1*<jS.1  V  j»»> 

X1 

Now  we  have  to  l  nd  transformations  arg  and  arg  which  enable  (195) 
to  equal  die  term  in  brackets  in  (194).  Thus  we  equate 
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2  M-l 


(  S  w.Xj-c )  *  (2  w4x4-arg, ) 
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and 


N  M  , 

S  ^  NiVV 

i«l  jsl  J  x 


N  M-l  2 

S  <  S  a«V*“g1  V 
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(197) 


From  (196) 


arg2 


C"WM*M 


(190) 


If  we  fix  i  in  (197)  we  obtain 

”*1  =  b-*M*(M)  I1”* 

where  s  (*jMt  a2M'****aNM^  aa^  **  transpose  of  the  M**1 
column  in  the  A  matrix.  Using  (198)  and  (199)  we  can  write  the  recurrence 
relation  (192). 

Computer  Programming 

The  computational  results  of  this  paper  were  obtained  using  the  IBM 
1620  and  IBM  7094  digital  computers*  The  1620  was  available  at  Utah 
State  University.  The  7094  was  available  through  the  Western  Data 
Processing  Center  at  UCLA  in  Los  Angeles,  California* 
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Two  facets  of  computer  programming  should  be  mentioned.  First, 
considerable  efforf  was  made  to  ensure  the  correct  estimation  of  image 
accuracy.  (This  effort  is  discussed  under  Image  and  A  Matrix  Accuracy 
in  the  Simulated  Object  Restorations  section. )  A  similar  effort  was 
made  to  ensure  that  the  solutions  were  correct.  A  good  check  of  these 
solutions  was  available  by  comparing  the  last  solution  vector  component 
*M  **  both  the  mAtri*  Aversion  solutions  and  the  corresponding  dynamic 
programming  solutions.  Second,  as  a  further  check  on  the  solution 
accuracy  and  to  study  the  inversion  process,  the  matrix  inverse  was 

checked  for  each  solution.  In  all  of  the  cases  shown  when  tr  2  was 

2  * 
finite,  <r^  was  assumed  small  enough  »  that  no  special  programming 

sophistication  was  required  to  ensure  a  correct  matrix  inverse.  When 

°x  Wa*  °*  very  large,  results  were  obtained  on  the  1620  which 

allowed  the  use  of  excr  isive  accuracy  (up  to  28  significant  digits).  This 

e.  .cessive  accuracy  enabled  results  to  be  obtained,  as  previously  presented, 

when  tiie  A  matrix  was  nearly  singular. 

A®  811  example  of  the  computer  programming  necessary  to  simulate 

the  restoration  of  optical  objects,  we  present  the  following  computer 

« 

program  which  represents  roughly  50  percent  of  the  computational  effort. 
This  computer  program  was  used  to  simulate  the  matrix  solutions  in  the 
restoration  problem  for  the  sources  composed  of  cosine  pulses  from 
cos  (kir(or-<p)]  as  shown  in  (166)  and  068).  The  random  number  subroutine 


IWIWWiiWIWWinnillllBilwiili^  iirrhiiiiiiMWieiiiw»wwe«MW»rt*wiwweiwweeMewwweewwiwp>w>iweeewweww«wgw^^ 


203 


was  obtained  from  the  Western  Data  Processing  Center  in  Los  Angeles* 
The  symbols  used  and  the  computational  setjuoD  e  are  indicated  as  the 
prog  ’am  proceeds. 
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